Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_SerialDenseSVD.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42*/
43
46
47//=============================================================================
49 : Epetra_Object("Epetra::SerialDenseSVD"),
53 Transpose_(false),
54 Factored_(false),
55 Solved_(false),
56 Inverted_(false),
57 TRANS_('N'),
58 M_(0),
59 N_(0),
60 Min_MN_(0),
61 NRHS_(0),
62 LDA_(0),
63 LDAI_(0),
64 LDB_(0),
65 LDX_(0),
66 INFO_(0),
67 LWORK_(0),
68 IWORK_(0),
69 ANORM_(0.0),
70 Matrix_(0),
71 LHS_(0),
72 RHS_(0),
73 Inverse_(0),
74 A_(0),
75 AI_(0),
76 WORK_(0),
77 U_(0),
78 S_(0),
79 Vt_(0),
80 B_(0),
81 X_(0),
82 UseTranspose_(false)
83{
87}
88//=============================================================================
93//=============================================================================
95{
96 IWORK_ = 0;
97// FERR_ = 0;
98// BERR_ = 0;
99// Factor_ =0;
100 Inverse_ =0;
101// AF_ = 0;
102 AI_ = 0;
103// IPIV_ = 0;
104 WORK_ = 0;
105// R_ = 0;
106// C_ = 0;
107 U_ = 0;
108 S_ = 0;
109 Vt_ = 0;
110 INFO_ = 0;
111 LWORK_ = 0;
112}
113//=============================================================================
115{
116 if (U_ !=0) {delete[] U_; U_ = 0;}
117 if (S_ !=0) {delete[] S_; S_ = 0;}
118 if (Vt_ !=0) {delete[] Vt_; Vt_ = 0;}
119 if (IWORK_ != 0) {delete [] IWORK_; IWORK_ = 0;}
120// if (FERR_ != 0) {delete [] FERR_; FERR_ = 0;}
121// if (BERR_ != 0) {delete [] BERR_; BERR_ = 0;}
122// if (Factor_ != Matrix_ && Factor_ != 0) {delete Factor_; Factor_ = 0;}
123// if (Factor_ !=0) Factor_ = 0;
124 if (Inverse_ != 0) {delete Inverse_; Inverse_ = 0;}
125// if (AF_ !=0) AF_ = 0;
126 if (AI_ !=0) AI_ = 0;
127// if (IPIV_ != 0) {delete [] IPIV_;IPIV_ = 0;}
128 if (WORK_ != 0) {delete [] WORK_;WORK_ = 0;}
129// if (R_ != 0 && R_ != C_) {delete [] R_;R_ = 0;}
130// if (R_ != 0) R_ = 0;
131// if (C_ != 0) {delete [] C_;C_ = 0;}
132 INFO_ = 0;
133 LWORK_ = 0;
134}
135//=============================================================================
137{
138 DeleteArrays();
139 ResetVectors();
140 Matrix_ = 0;
141 Inverse_ = 0;
142// Factor_ = 0;
143// A_Equilibrated_ = false;
144 Factored_ = false;
145 Inverted_ = false;
146 M_ = 0;
147 N_ = 0;
148 Min_MN_ = 0;
149 LDA_ = 0;
150// LDAF_ = 0;
151 LDAI_ = 0;
152 ANORM_ = -1.0;
153// RCOND_ = -1.0;
154// ROWCND_ = -1.0;
155// COLCND_ = -1.0;
156// AMAX_ = -1.0;
157 A_ = 0;
158
159 if( U_ ) { delete [] U_; U_ = 0; }
160 if( S_ ) { delete [] S_; S_ = 0; }
161 if( Vt_ ) { delete [] Vt_; Vt_ = 0; }
162}
163//=============================================================================
165 ResetMatrix();
166 Matrix_ = &A_in;
167// Factor_ = &A_in;
168 M_ = A_in.M();
169 N_ = A_in.N();
171 LDA_ = A_in.LDA();
172// LDAF_ = LDA_;
173 A_ = A_in.A();
174// AF_ = A_in.A();
175 return(0);
176}
177//=============================================================================
179{
180 LHS_ = 0;
181 RHS_ = 0;
182 B_ = 0;
183 X_ = 0;
184// ReciprocalConditionEstimated_ = false;
185// SolutionRefined_ = false;
186 Solved_ = false;
187// SolutionErrorsEstimated_ = false;
188// B_Equilibrated_ = false;
189 NRHS_ = 0;
190 LDB_ = 0;
191 LDX_ = 0;
192}
193//=============================================================================
195{
196 if (B_in.M()!=X_in.M() || B_in.N() != X_in.N()) EPETRA_CHK_ERR(-1);
197 if (B_in.A()==0) EPETRA_CHK_ERR(-2);
198 if (B_in.LDA()<1) EPETRA_CHK_ERR(-3);
199 if (X_in.A()==0) EPETRA_CHK_ERR(-4);
200 if (X_in.LDA()<1) EPETRA_CHK_ERR(-5);
201
202 ResetVectors();
203 LHS_ = &X_in;
204 RHS_ = &B_in;
205 NRHS_ = B_in.N();
206
207 B_ = B_in.A();
208 LDB_ = B_in.LDA();
209 X_ = X_in.A();
210 LDX_ = X_in.LDA();
211 return(0);
212}
213//=============================================================================
215
216 int ierr = 0;
217
218 ANORM_ = Matrix_->OneNorm(); // Compute 1-Norm of A
219
220 //allocate U_, S_, and Vt_ if not already done
221 if(U_==0)
222 {
223 U_ = new double[M_*N_];
224 S_ = new double[M_];
225 Vt_ = new double[M_*N_];
226 }
227 else //zero them out
228 {
229 for( int i = 0; i < M_; ++i ) S_[i]=0.0;
230 for( int i = 0; i < M_*N_; ++i )
231 {
232 U_[i]=0.0;
233 Vt_[i]=0.0;
234 }
235 }
236
237// if (Equilibrate_) ierr = EquilibrateMatrix();
238
239 if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
240
241 //allocate temp work space
242 int lwork = 5*M_;
243 double *work = new double[lwork];
244 char job = 'A';
245
246 //create temporary matrix to avoid writeover of original
248 GESVD( job, job, M_, N_, tempMat.A(), LDA_, S_, U_, N_, Vt_, M_, work, &lwork, &INFO_ );
249
250 delete [] work;
251
252 Factored_ = true;
253 double DN = N_;
254 UpdateFlops(2.0*(DN*DN*DN)/3.0);
255
257 return(0);
258}
259
260//=============================================================================
262
263 //FOR NOW, ONLY ALLOW SOLVES IF INVERTED!!!!
264 //NO REFINEMENT!!!
265 //NO EQUILIBRATION!!!
266
267 // We will call one of four routines depending on what services the user wants and
268 // whether or not the matrix has been inverted or factored already.
269 //
270 // If the matrix has been inverted, use DGEMM to compute solution.
271 // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
272 // call the X interface.
273 // Otherwise, if the matrix is already factored we will call the TRS interface.
274 // Otherwise, if the matrix is unfactored we will call the SV interface.
275
276
277/*
278 if (Equilibrate_) {
279 ierr = EquilibrateRHS();
280 B_Equilibrated_ = true;
281 }
282 EPETRA_CHK_ERR(ierr);
283 if (A_Equilibrated_ && !B_Equilibrated_) EPETRA_CHK_ERR(-1); // Matrix and vectors must be similarly scaled
284 if (!A_Equilibrated_ && B_Equilibrated_) EPETRA_CHK_ERR(-2);
285 if (B_==0) EPETRA_CHK_ERR(-3); // No B
286 if (X_==0) EPETRA_CHK_ERR(-4); // No B
287
288 if (ShouldEquilibrate() && !A_Equilibrated_) ierr = 1; // Warn that the system should be equilibrated.
289*/
290
291 double DN = N_;
292 double DNRHS = NRHS_;
293 if (Inverted()) {
294
295 if (B_==X_) EPETRA_CHK_ERR(-100); // B and X must be different for this case
296
297 GEMM(TRANS_, 'N', N_, NRHS_, N_, 1.0, AI_, LDAI_, B_, LDB_, 0.0, X_, LDX_);
298 if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
299 UpdateFlops(2.0*DN*DN*DNRHS);
300 Solved_ = true;
301 }
302 else EPETRA_CHK_ERR(-101); //Currently, for solve must have inverse
303/*
304 else {
305
306 if (!Factored()) Factor(); // Matrix must be factored
307
308 if (B_!=X_) *LHS_ = *RHS_; // Copy B to X if needed
309 GETRS(TRANS_, N_, NRHS_, AF_, LDAF_, IPIV_, X_, LDX_, &INFO_);
310 if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
311 UpdateFlops(2.0*DN*DN*DNRHS);
312 Solved_ = true;
313
314 }
315
316 int ierr1=0;
317 if (RefineSolution_ && !Inverted()) ierr1 = ApplyRefinement();
318 if (ierr1!=0) EPETRA_CHK_ERR(ierr1)
319 else
320 EPETRA_CHK_ERR(ierr);
321
322 if (Equilibrate_) ierr1 = UnequilibrateLHS();
323 EPETRA_CHK_ERR(ierr1);
324*/
325 return(0);
326}
327/*
328//=============================================================================
329int Epetra_SerialDenseSVD::ApplyRefinement(void)
330{
331 double DN = N_;
332 double DNRHS = NRHS_;
333 if (!Solved()) EPETRA_CHK_ERR(-100); // Must have an existing solution
334 if (A_==AF_) EPETRA_CHK_ERR(-101); // Cannot apply refine if no original copy of A.
335
336 if (FERR_ != 0) delete [] FERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
337 FERR_ = new double[NRHS_];
338 if (BERR_ != 0) delete [] BERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
339 BERR_ = new double[NRHS_];
340 AllocateWORK();
341 AllocateIWORK();
342
343 GERFS(TRANS_, N_, NRHS_, A_, LDA_, AF_, LDAF_, IPIV_,
344 B_, LDB_, X_, LDX_, FERR_, BERR_,
345 WORK_, IWORK_, &INFO_);
346
347
348 SolutionErrorsEstimated_ = true;
349 ReciprocalConditionEstimated_ = true;
350 SolutionRefined_ = true;
351
352 UpdateFlops(2.0*DN*DN*DNRHS); // Not sure of count
353
354 EPETRA_CHK_ERR(INFO_);
355 return(0);
356
357}
358
359//=============================================================================
360int Epetra_SerialDenseSVD::ComputeEquilibrateScaling(void) {
361 if (R_!=0) return(0); // Already computed
362
363 double DM = M_;
364 double DN = N_;
365 R_ = new double[M_];
366 C_ = new double[N_];
367
368 GEEQU (M_, N_, AF_, LDAF_, R_, C_, &ROWCND_, &COLCND_, &AMAX_, &INFO_);
369 if (INFO_ != 0) EPETRA_CHK_ERR(INFO_);
370
371 if (COLCND_<0.1 || ROWCND_<0.1 || AMAX_ < Epetra_Underflow || AMAX_ > Epetra_Overflow) ShouldEquilibrate_ = true;
372
373 UpdateFlops(4.0*DM*DN);
374
375 return(0);
376}
377
378//=============================================================================
379int Epetra_SerialDenseSVD::EquilibrateMatrix(void)
380{
381 int i, j;
382 int ierr = 0;
383
384 double DN = N_;
385 double DM = M_;
386
387 if (A_Equilibrated_) return(0); // Already done
388 if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute R and C if needed
389 if (ierr!=0) EPETRA_CHK_ERR(ierr);
390 if (A_==AF_) {
391 double * ptr;
392 for (j=0; j<N_; j++) {
393 ptr = A_ + j*LDA_;
394 double s1 = C_[j];
395 for (i=0; i<M_; i++) {
396 *ptr = *ptr*s1*R_[i];
397 ptr++;
398 }
399 }
400 UpdateFlops(2.0*DM*DN);
401 }
402 else {
403 double * ptr;
404 double * ptr1;
405 for (j=0; j<N_; j++) {
406 ptr = A_ + j*LDA_;
407 ptr1 = AF_ + j*LDAF_;
408 double s1 = C_[j];
409 for (i=0; i<M_; i++) {
410 *ptr = *ptr*s1*R_[i];
411 ptr++;
412 *ptr1 = *ptr1*s1*R_[i];
413 ptr1++;
414 }
415 }
416 UpdateFlops(4.0*DM*DN);
417 }
418
419 A_Equilibrated_ = true;
420
421 return(0);
422}
423
424//=============================================================================
425int Epetra_SerialDenseSVD::EquilibrateRHS(void)
426{
427 int i, j;
428 int ierr = 0;
429
430 if (B_Equilibrated_) return(0); // Already done
431 if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute R and C if needed
432 if (ierr!=0) EPETRA_CHK_ERR(ierr);
433
434 double * R = R_;
435 if (Transpose_) R = C_;
436
437 double * ptr;
438 for (j=0; j<NRHS_; j++) {
439 ptr = B_ + j*LDB_;
440 for (i=0; i<M_; i++) {
441 *ptr = *ptr*R[i];
442 ptr++;
443 }
444 }
445
446
447 B_Equilibrated_ = true;
448 UpdateFlops((double) N_*(double) NRHS_);
449
450 return(0);
451}
452
453//=============================================================================
454int Epetra_SerialDenseSVD::UnequilibrateLHS(void)
455{
456 int i, j;
457
458 if (!B_Equilibrated_) return(0); // Nothing to do
459
460 double * C = C_;
461 if (Transpose_) C = R_;
462
463 double * ptr;
464 for (j=0; j<NRHS_; j++) {
465 ptr = X_ + j*LDX_;
466 for (i=0; i<N_; i++) {
467 *ptr = *ptr*C[i];
468 ptr++;
469 }
470 }
471
472
473 UpdateFlops((double) N_ *(double) NRHS_);
474
475 return(0);
476}
477*/
478
479//=============================================================================
480int Epetra_SerialDenseSVD::Invert( double rthresh, double athresh )
481{
482 if (!Factored()) Factor(); // Need matrix factored.
483
484 //apply threshold
485 double thresh = S_[0]*rthresh + athresh;
486 int num_replaced = 0;
487 for( int i = 0; i < M_; ++i )
488 if( S_[i] < thresh )
489 {
490//cout << num_replaced << thresh << " " << S_[0] << " " << S_[i] << std::endl;
491// S_[i] = thresh;
492 S_[i] = 0.0;
493 ++num_replaced;
494 }
495
496 //scale cols of U_ with reciprocal singular values
497 double *p = U_;
498 for( int i = 0; i < N_; ++i )
499 {
500 double scale = 0.0;
501 if( S_[i] ) scale = 1./S_[i];
502 for( int j = 0; j < M_; ++j ) *p++ *= scale;
503 }
504
505 //create new Inverse_ if necessary
506 if( Inverse_ == 0 )
507 {
509 Inverse_->Shape( N_, M_ );
510 AI_ = Inverse_->A();
511 LDAI_ = Inverse_->LDA();
512 }
513/*
514 else //zero it out
515 {
516 for( int i = 0; i < Inverse_->M(); ++i )
517 for( int j = 0; j < Inverse_->N(); ++j )
518 (*Inverse_)(i,j) = 0.0;
519 }
520*/
521
522 GEMM( 'T', 'T', M_, M_, M_, 1.0, Vt_, M_, U_, M_, 0.0, AI_, M_ );
523
524 double DN = N_;
525 UpdateFlops((DN*DN*DN));
526 Inverted_ = true;
527 Factored_ = false;
528
530 return(num_replaced);
531}
532
533/*
534//=============================================================================
535int Epetra_SerialDenseSVD::ReciprocalConditionEstimate(double & Value)
536{
537 int ierr = 0;
538 if (ReciprocalConditionEstimated()) {
539 Value = RCOND_;
540 return(0); // Already computed, just return it.
541 }
542
543 if (ANORM_<0.0) ANORM_ = Matrix_->OneNorm();
544 if (!Factored()) ierr = Factor(); // Need matrix factored.
545 if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
546
547 AllocateWORK();
548 AllocateIWORK();
549 // We will assume a one-norm condition number
550 GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, WORK_, IWORK_, &INFO_);
551 ReciprocalConditionEstimated_ = true;
552 Value = RCOND_;
553 UpdateFlops(2*N_*N_); // Not sure of count
554 EPETRA_CHK_ERR(INFO_);
555 return(0);
556}
557*/
558
559//=============================================================================
560void Epetra_SerialDenseSVD::Print(std::ostream& os) const {
561
562 if (Matrix_!=0) os << *Matrix_;
563// if (Factor_!=0) os << *Factor_;
564 if (S_!=0) for( int i = 0; i < M_; ++i ) std::cout << "(" << i << "," << S_[i] << ")\n";
565 if (Inverse_!=0) os << *Inverse_;
566 if (LHS_!=0) os << *LHS_;
567 if (RHS_!=0) os << *RHS_;
568
569}
#define EPETRA_MIN(x, y)
#define EPETRA_CHK_ERR(a)
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition Epetra_BLAS.h:70
void GEMM(const char TRANSA, const char TRANSB, const int M, const int N, const int K, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
Epetra_BLAS matrix-matrix multiply function (SGEMM)
Epetra_CompObject: Functionality and data that is common to all computational classes.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
Epetra_LAPACK: The Epetra LAPACK Wrapper Class.
void GESVD(const char JOBU, const char JOBVT, const int M, const int N, float *A, const int LDA, float *S, float *U, const int LDU, float *VT, const int LDVT, float *WORK, const int *LWORK, int *INFO) const
Epetra_LAPACK wrapper for computing the singular value decomposition (SGESVD)
Epetra_Object: The base Epetra class.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
double * A() const
Returns pointer to the this matrix.
int LDA() const
Returns the leading dimension of the this matrix.
virtual double OneNorm() const
Computes the 1-Norm of the this matrix (identical to NormOne() method).
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
int M() const
Returns row dimension of system.
int N() const
Returns column dimension of system.
Epetra_SerialDenseMatrix * RHS_
Epetra_SerialDenseMatrix * Inverse_
Epetra_SerialDenseMatrix * Matrix_
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
virtual int Invert(double rthresh=0.0, double athresh=0.0)
Inverts the this matrix.
virtual ~Epetra_SerialDenseSVD()
Epetra_SerialDenseSVD destructor.
Epetra_SerialDenseMatrix * LHS_
int SetMatrix(Epetra_SerialDenseMatrix &A)
Sets the pointers for coefficient matrix.
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
Epetra_SerialDenseSVD()
Default constructor; matrix should be set using SetMatrix(), LHS and RHS set with SetVectors().
virtual void Print(std::ostream &os) const
Print service methods; defines behavior of ostream << operator.
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).