00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "Epetra_config.h"
00043 #include "LOCA_Epetra_LowRankUpdateRowMatrix.H"
00044
00045 #include "Epetra_Vector.h"
00046 #include "Epetra_Map.h"
00047 #include "Epetra_Comm.h"
00048 #include "LOCA_GlobalData.H"
00049 #include "LOCA_ErrorCheck.H"
00050
00051 LOCA::Epetra::LowRankUpdateRowMatrix::
00052 LowRankUpdateRowMatrix(
00053 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00054 const Teuchos::RCP<Epetra_RowMatrix>& jacRowMatrix,
00055 const Teuchos::RCP<Epetra_MultiVector>& U_multiVec,
00056 const Teuchos::RCP<Epetra_MultiVector>& V_multiVec,
00057 bool setup_for_solve,
00058 bool include_UV_terms) :
00059 LOCA::Epetra::LowRankUpdateOp(global_data, jacRowMatrix, U_multiVec,
00060 V_multiVec, setup_for_solve),
00061 J_rowMatrix(jacRowMatrix),
00062 nonconst_U(U_multiVec),
00063 nonconst_V(V_multiVec),
00064 includeUV(include_UV_terms),
00065 m(U_multiVec->NumVectors()),
00066 U_map(U_multiVec->Map()),
00067 V_map(V_multiVec->Map()),
00068 row_map(jacRowMatrix->RowMatrixRowMap())
00069 {
00070 }
00071
00072 LOCA::Epetra::LowRankUpdateRowMatrix::
00073 ~LowRankUpdateRowMatrix()
00074 {
00075 }
00076
00077 const Epetra_BlockMap &
00078 LOCA::Epetra::LowRankUpdateRowMatrix::
00079 Map() const
00080 {
00081 return J_rowMatrix->Map();
00082 }
00083
00084 int
00085 LOCA::Epetra::LowRankUpdateRowMatrix::
00086 NumMyRowEntries(int MyRow, int & NumEntries) const
00087 {
00088 return J_rowMatrix->NumMyRowEntries(MyRow, NumEntries);
00089 }
00090
00091 int
00092 LOCA::Epetra::LowRankUpdateRowMatrix::
00093 MaxNumEntries() const
00094 {
00095 return J_rowMatrix->MaxNumEntries();
00096 }
00097
00098 int
00099 LOCA::Epetra::LowRankUpdateRowMatrix::
00100 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values,
00101 int * Indices) const
00102 {
00103
00104 int res = J_rowMatrix->ExtractMyRowCopy(MyRow, Length, NumEntries, Values,
00105 Indices);
00106
00107 if (!includeUV)
00108 return res;
00109
00110
00111 int jac_row_gid = row_map.GID(MyRow);
00112 int u_row_lid = U_map.LID(jac_row_gid);
00113 for (int col=0; col<NumEntries; col++) {
00114 int jac_col_gid = Indices[col];
00115 if (V_map.MyGID(jac_col_gid)) {
00116 int v_row_lid = V_map.LID(jac_col_gid);
00117 Values[col] += computeUV(u_row_lid, v_row_lid);
00118 }
00119 }
00120
00121 return res;
00122 }
00123
00124 int
00125 LOCA::Epetra::LowRankUpdateRowMatrix::
00126 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00127 {
00128
00129 int res = J_rowMatrix->ExtractDiagonalCopy(Diagonal);
00130
00131 if (!includeUV)
00132 return res;
00133
00134
00135 int numMyRows = J_rowMatrix->NumMyRows();
00136 for (int row=0; row<numMyRows; row++) {
00137 int jac_row_gid = row_map.GID(row);
00138 int u_row_lid = U_map.LID(jac_row_gid);
00139 int v_row_lid = V_map.LID(jac_row_gid);
00140 Diagonal[row] += computeUV(u_row_lid, v_row_lid);
00141 }
00142
00143 return res;
00144 }
00145
00146 int
00147 LOCA::Epetra::LowRankUpdateRowMatrix::
00148 Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00149 {
00150
00151 int res = J_rowMatrix->Multiply(TransA, X, Y);
00152
00153
00154 int numVecs = X.NumVectors();
00155
00156
00157 if (tmpMat == Teuchos::null || tmpMat->NumVectors() != numVecs) {
00158 tmpMat = Teuchos::rcp(new Epetra_MultiVector(localMap, numVecs, false));
00159 }
00160
00161 if (!TransA) {
00162
00163
00164 tmpMat->Multiply('T', 'N', 1.0, *V, X, 0.0);
00165
00166
00167 Y.Multiply('N', 'N', 1.0, *U, *tmpMat, 1.0);
00168
00169 }
00170 else {
00171
00172
00173 tmpMat->Multiply('T', 'N', 1.0, *U, X, 0.0);
00174
00175
00176 Y.Multiply('N', 'N', 1.0, *V, *tmpMat, 1.0);
00177
00178 }
00179
00180 return res;
00181 }
00182
00183 int
00184 LOCA::Epetra::LowRankUpdateRowMatrix::
00185 Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X,
00186 Epetra_MultiVector& Y) const
00187 {
00188
00189 return J_rowMatrix->Solve(Upper, Trans, UnitDiagonal, X, Y);
00190 }
00191
00192 int
00193 LOCA::Epetra::LowRankUpdateRowMatrix::
00194 InvRowSums(Epetra_Vector& x) const
00195 {
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218 int res;
00219
00220 if (!includeUV)
00221 res = J_rowMatrix->InvRowSums(x);
00222 else {
00223 Epetra_Vector ones(V_map,false);
00224 ones.PutScalar(1.0);
00225 res = Multiply(true, ones, x);
00226 x.Reciprocal(x);
00227 }
00228
00229 return res;
00230 }
00231
00232 int
00233 LOCA::Epetra::LowRankUpdateRowMatrix::
00234 LeftScale(const Epetra_Vector& x)
00235 {
00236
00237
00238 int res;
00239 if (J_rowMatrix->UseTranspose())
00240 res = J_rowMatrix->RightScale(x);
00241 else
00242 res = J_rowMatrix->LeftScale(x);
00243
00244
00245 for (int j=0; j<m; j++)
00246 (*nonconst_U)(j)->Multiply(1.0, x, *((*nonconst_U)(j)), 0.0);
00247
00248 return res;
00249 }
00250
00251 int
00252 LOCA::Epetra::LowRankUpdateRowMatrix::
00253 InvColSums(Epetra_Vector& x) const
00254 {
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 int res;
00278
00279 if (!includeUV)
00280 res = J_rowMatrix->InvColSums(x);
00281 else {
00282 Epetra_Vector ones(V_map,false);
00283 ones.PutScalar(1.0);
00284 res = Multiply(false, ones, x);
00285 x.Reciprocal(x);
00286 }
00287
00288 return res;
00289 }
00290
00291 int
00292 LOCA::Epetra::LowRankUpdateRowMatrix::
00293 RightScale(const Epetra_Vector& x)
00294 {
00295
00296 int res;
00297 if (J_rowMatrix->UseTranspose())
00298 res = J_rowMatrix->LeftScale(x);
00299 else
00300 res = J_rowMatrix->RightScale(x);
00301
00302
00303 for (int j=0; j<m; j++)
00304 (*nonconst_V)(j)->Multiply(1.0, x, *((*nonconst_V)(j)), 0.0);
00305
00306 return res;
00307 }
00308
00309 bool
00310 LOCA::Epetra::LowRankUpdateRowMatrix::
00311 Filled() const
00312 {
00313 return J_rowMatrix->Filled();
00314 }
00315
00316 double
00317 LOCA::Epetra::LowRankUpdateRowMatrix::
00318 NormInf() const
00319 {
00320 if (!includeUV)
00321 return J_rowMatrix->NormInf();
00322
00323 Epetra_Vector tmp(row_map);
00324 InvRowSums(tmp);
00325 tmp.Reciprocal(tmp);
00326
00327 double val;
00328 tmp.MaxValue(&val);
00329
00330 return val;
00331 }
00332
00333 double
00334 LOCA::Epetra::LowRankUpdateRowMatrix::
00335 NormOne() const
00336 {
00337 if (!includeUV)
00338 return J_rowMatrix->NormOne();
00339
00340 Epetra_Vector tmp(row_map);
00341 InvColSums(tmp);
00342 tmp.Reciprocal(tmp);
00343
00344 double val;
00345 tmp.MaxValue(&val);
00346
00347 return val;
00348 }
00349
00350 int
00351 LOCA::Epetra::LowRankUpdateRowMatrix::
00352 NumGlobalNonzeros() const
00353 {
00354 return J_rowMatrix->NumGlobalNonzeros();
00355 }
00356
00357 int
00358 LOCA::Epetra::LowRankUpdateRowMatrix::
00359 NumGlobalRows() const
00360 {
00361 return J_rowMatrix->NumGlobalRows();
00362 }
00363
00364 int
00365 LOCA::Epetra::LowRankUpdateRowMatrix::
00366 NumGlobalCols() const
00367 {
00368 return J_rowMatrix->NumGlobalCols();
00369 }
00370
00371 int
00372 LOCA::Epetra::LowRankUpdateRowMatrix::
00373 NumGlobalDiagonals() const
00374 {
00375 return J_rowMatrix->NumGlobalDiagonals();
00376 }
00377
00378 int
00379 LOCA::Epetra::LowRankUpdateRowMatrix::
00380 NumMyNonzeros() const
00381 {
00382 return J_rowMatrix->NumMyNonzeros();
00383 }
00384
00385 int
00386 LOCA::Epetra::LowRankUpdateRowMatrix::
00387 NumMyRows() const
00388 {
00389 return J_rowMatrix->NumMyRows();
00390 }
00391
00392 int
00393 LOCA::Epetra::LowRankUpdateRowMatrix::
00394 NumMyCols() const
00395 {
00396 return J_rowMatrix->NumMyCols();
00397 }
00398
00399 int
00400 LOCA::Epetra::LowRankUpdateRowMatrix::
00401 NumMyDiagonals() const
00402 {
00403 return J_rowMatrix->NumMyDiagonals();
00404 }
00405
00406 bool
00407 LOCA::Epetra::LowRankUpdateRowMatrix::
00408 LowerTriangular() const
00409 {
00410 return J_rowMatrix->LowerTriangular();
00411 }
00412
00413 bool
00414 LOCA::Epetra::LowRankUpdateRowMatrix::
00415 UpperTriangular() const
00416 {
00417 return J_rowMatrix->UpperTriangular();
00418 }
00419
00420 const Epetra_Map &
00421 LOCA::Epetra::LowRankUpdateRowMatrix::
00422 RowMatrixRowMap() const
00423 {
00424 return J_rowMatrix->RowMatrixRowMap();
00425 }
00426
00427 const Epetra_Map &
00428 LOCA::Epetra::LowRankUpdateRowMatrix::
00429 RowMatrixColMap() const
00430 {
00431 return J_rowMatrix->RowMatrixColMap();
00432 }
00433
00434 const Epetra_Import *
00435 LOCA::Epetra::LowRankUpdateRowMatrix::
00436 RowMatrixImporter() const
00437 {
00438 return J_rowMatrix->RowMatrixImporter();
00439 }
00440
00441 double
00442 LOCA::Epetra::LowRankUpdateRowMatrix::
00443 computeUV(int u_row_lid, int v_row_lid) const
00444 {
00445 double val = 0.0;
00446
00447
00448 for (int k=0; k<m; k++)
00449 val += (*U)[k][u_row_lid] * (*V)[k][v_row_lid];
00450
00451 return val;
00452 }