RBGen  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
RBGen_ISVDSingle.cpp
1 #include "RBGen_ISVDSingle.h"
3 #include "Epetra_Comm.h"
4 
5 namespace RBGen {
6 
8 
10  // free pointer to original data set
11  // switch it to new data set, momentarily
12  A_ = update_ss;
13 
14  // we may be calling update basis from the beginning
15  // we hope that V_ (if it exists) is tall enough, we will check this
16 
17  // starting with the current factorization, make a single pass
18  const int oldNumCols = numProc_;
19  const int newNumCols = A_->NumVectors();
20  TEUCHOS_TEST_FOR_EXCEPTION(oldNumCols+newNumCols > maxNumCols_, std::invalid_argument,
21  "RBGen::ISVSingle::updateBasis(): number of snapshots to process has exceeded specified maximum humber of columns.");
22  while (numProc_ < oldNumCols+newNumCols) {
23  // determine lup
24  int lup;
25  if (curRank_ == 0) {
26  // first step
27  lup = startRank_;
28  }
29  else {
30  // this value minimizes overall complexity for a UDV factorization, assuming fixed rank
31  lup = (int)(curRank_ / Teuchos::ScalarTraits<double>::squareroot(2.0));
32  }
33  // now cap lup via lmin,lmax,maxBasisSize
34  // want lup >= lmin
35  // need lup <= numCols - numProc
36  // lup <= lmax
37  // lup <= maxBasisSize - curRank
38  lup = (lup < lmin_ ? lmin_ : lup);
39  lup = (lup > oldNumCols+newNumCols - numProc_ ? oldNumCols+newNumCols - numProc_ : lup);
40  lup = (lup > lmax_ ? lmax_ : lup);
41  lup = (lup > maxBasisSize_ - curRank_ ? maxBasisSize_ - curRank_ : lup);
42 
43  // get view of new vectors
46  Aplus = Teuchos::rcp( new Epetra_MultiVector(::View,*A_,numProc_-oldNumCols,lup));
47  Unew = Teuchos::rcp( new Epetra_MultiVector(::View,*U_,curRank_,lup));
48  // put them in U
49  *Unew = *Aplus;
50  // clear the views
51  Unew = Teuchos::null;
52  Aplus = Teuchos::null;
53 
54  // perform the incremental step
55  incStep(lup);
56  }
57 
58  //
59  // we can't compute residuals, because we don't have old and new snapshots
60  for (int i=0; i<curRank_; i++) {
61  resNorms_[i] = 0;
62  }
63 
64  // print out some info
65  const Epetra_Comm *comm = &A_->Comm();
66  if (comm->MyPID() == 0 && verbLevel_ >= 1) {
67  std::cout
68  << "------------- ISVDSingle::updateBasis() -----------" << std::endl
69  << "| Current rank: " << curRank_ << std::endl
70  << "| Current sigmas: " << std::endl;
71  for (int i=0; i<curRank_; i++) {
72  std::cout << "| " << sigma_[i] << std::endl;
73  }
74  }
75 
76  return;
77  }
78 
79  void ISVDSingle::makePass() {
80  // ISVDSingle only makes a single pass
81  TEUCHOS_TEST_FOR_EXCEPTION(maxNumPasses_ != 1,std::logic_error,
82  "RBGen::ISVDSingle::makePass(): Max Num Passes should be 1, but is not.");
83  // did we already make our one pass? we can't afford to run this again
84  if (curNumPasses_ > 0) return;
85  const int numCols = A_->NumVectors();
86  while (numProc_ < numCols) {
87  // determine lup
88  int lup;
89  if (curRank_ == 0) {
90  // first step
91  lup = startRank_;
92  }
93  else {
94  // this value minimizes overall complexity, assuming fixed rank
95  lup = (int)(curRank_ / Teuchos::ScalarTraits<double>::squareroot(2.0));
96  }
97  // now cap lup via lmin,lmax,maxBasisSize
98  // want lup >= lmin
99  // need lup <= numCols - numProc
100  // lup <= lmax
101  // lup <= maxBasisSize - curRank
102  lup = (lup < lmin_ ? lmin_ : lup);
103  lup = (lup > numCols - numProc_ ? numCols - numProc_ : lup);
104  lup = (lup > lmax_ ? lmax_ : lup);
105  lup = (lup > maxBasisSize_ - curRank_ ? maxBasisSize_ - curRank_ : lup);
106 
107  // get view of new vectors
110  Aplus = Teuchos::rcp( new Epetra_MultiVector(::View,*A_,numProc_,lup));
111  Unew = Teuchos::rcp( new Epetra_MultiVector(::View,*U_,curRank_,lup));
112  // put them in U
113  *Unew = *Aplus;
114  Unew = Teuchos::null;
115  Aplus = Teuchos::null;
116 
117  // perform the incremental step
118  incStep(lup);
119  }
120 
121  //
122  // compute the new residuals
123  // we know that A V = U S
124  // if, in addition, A^T U = V S, then have singular subspaces
125  // check residuals A^T U - V S, scaling the i-th column by sigma[i]
126  //
127  {
128  Epetra_LocalMap lclmap(A_->NumVectors(),0,A_->Comm());
129  Epetra_MultiVector ATU(lclmap,maxBasisSize_,false);
130 
131  // we know that A V = U S
132  // if, in addition, A^T U = V S, then have singular subspaces
133  // check residuals A^T U - V S, scaling the i-th column by sigma[i]
134  Epetra_MultiVector ATUlcl(::View,ATU,0,curRank_);
135  Epetra_MultiVector Ulcl(::View,*U_,0,curRank_);
136  Epetra_MultiVector Vlcl(::View,*V_,0,curRank_);
137  // compute A^T U
138  int info = ATUlcl.Multiply('T','N',1.0,*A_,Ulcl,0.0);
139  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
140  "RBGen::ISVDMultiCD::makePass(): Error calling Epetra_MultiVector::Multiply for A^T U.");
141  Epetra_LocalMap rankmap(curRank_,0,A_->Comm());
142  Epetra_MultiVector S(rankmap,curRank_,true);
143  for (int i=0; i<curRank_; i++) {
144  S[i][i] = sigma_[i];
145  }
146  // subtract V S from A^T U
147  info = ATUlcl.Multiply('N','N',-1.0,Vlcl,S,1.0);
148  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
149  "RBGen::ISVDMultiCD::computeBasis(): Error calling Epetra_MultiVector::Multiply for V S.");
150  resNorms_.resize(curRank_);
151  ATUlcl.Norm2(&resNorms_[0]);
152  // scale by sigmas
153  for (int i=0; i<curRank_; i++) {
154  if (sigma_[i] != 0.0) {
155  resNorms_[i] /= sigma_[i];
156  }
157  }
158 
159  }
160 
161 
162  // debugging checks
163  std::vector<double> errnorms(curRank_);
164  if (debug_) {
165  int info;
166  // Check that A V = U Sigma
167  // get pointers to current U and V, create workspace for A V - U Sigma
168  Epetra_MultiVector work(U_->Map(),curRank_,false),
169  curU(::View,*U_,0,curRank_),
170  curV(::View,*V_,0,curRank_);
171  // create local MV for sigmas
172  Epetra_LocalMap lclmap(curRank_,0,A_->Comm());
173  Epetra_MultiVector curS(lclmap,curRank_,true);
174  for (int i=0; i<curRank_; i++) {
175  curS[i][i] = sigma_[i];
176  }
177  info = work.Multiply('N','N',1.0,curU,curS,0.0);
178  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
179  "RBGen::ISVDMultiCD::makePass(): Error calling Epetra_MultiVector::Multiply() for debugging U S.");
180  info = work.Multiply('N','N',-1.0,*A_,curV,1.0);
181  TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
182  "RBGen::ISVDMultiCD::makePass(): Error calling Epetra_MultiVector::Multiply() for debugging U S - A V.");
183  work.Norm2(&errnorms[0]);
184  for (int i=0; i<curRank_; i++) {
185  if (sigma_[i] != 0.0) {
186  errnorms[i] /= sigma_[i];
187  }
188  }
189  }
190 
191 
192  // update pass counter
193  curNumPasses_++;
194 
195  // print out some info
196  const Epetra_Comm *comm = &A_->Comm();
197  if (comm->MyPID() == 0 && verbLevel_ >= 1) {
198  std::cout
199  << "------------- ISVDSingle::makePass() -----------" << std::endl
200  << "| Number of passes: " << curNumPasses_ << std::endl
201  << "| Current rank: " << curRank_ << std::endl
202  << "| Current sigmas: " << std::endl;
203  for (int i=0; i<curRank_; i++) {
204  std::cout << "| " << sigma_[i] << std::endl;
205  }
206  if (debug_) {
207  std::cout << "|DBG US-AV norms: " << std::endl;
208  for (int i=0; i<curRank_; i++) {
209  std::cout << "|DBG " << errnorms[i] << std::endl;
210  }
211  }
212  }
213 
214  return;
215  }
216 
221  )
222  {
223  maxNumPasses_ = 1;
224  }
225 
226 } // end of RBGen namespace
227 
228 
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
ISVDSingle()
Default constructor.
void updateBasis(const Teuchos::RCP< Epetra_MultiVector > &update_ss)
Update the current basis by appending new snapshots.
virtual int MyPID() const =0
void Initialize(const Teuchos::RCP< Teuchos::ParameterList > &params, const Teuchos::RCP< const Epetra_MultiVector > &init, const Teuchos::RCP< RBGen::FileIOHandler< Epetra_Operator > > &fileio=Teuchos::null)
Initialize the method with the given parameter list and snapshot set.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)