00001 #define max(a,b) ((a) >= (b) ? (a) : (b))
00002
00003 #include "myrlog.h"
00004
00005 #include "ml_include.h"
00006 #include "EpetraExt_MatrixMatrix.h"
00007 #include "Teuchos_ParameterList.hpp"
00008 #include "ml_epetra_preconditioner.h"
00009
00010 #include "OutputCapturer.h"
00011 #include "pbedefs.h"
00012
00013 #include "matrix_matrix.h"
00014 #include "directsolveroperator.h"
00015 #include "diagonalpreconditioner.h"
00016 #include "BlockPCGSolver.h"
00017
00018 #include "QuadraticEigsolvOperators.h"
00019
00020 QuadraticEigsolvOperators::QuadraticEigsolvOperators(FemaxMesh& mesh_mixed, Teuchos::ParameterList& params, const Epetra_Comm& MyComm) :
00021 Comm(MyComm),
00022 A(0),
00023 M(0),
00024 K(0),
00025 K11(0),
00026 K12(0),
00027 K22(0),
00028 H(0),
00029 H11(0),
00030 H12(0),
00031 H22(0),
00032 Y(0),
00033 Y11(0),
00034 C(0),
00035 Prec(0),
00036 h_solver_op_(0),
00037 HPrec(0),
00038 _sigma_K(0.0),
00039 params_(params),
00040 Aprec(0),
00041 AprecSolver(0),
00042 non(0),
00043 noe(0)
00044 {
00045 non = mesh_mixed.get_lagrange_mesh().get_num_global_mdofs(1);
00046 noe = mesh_mixed.get_nedelec_mesh().get_num_global_mdofs(1);
00047
00048 matrixAssembly(mesh_mixed);
00049 }
00050
00051 QuadraticEigsolvOperators::~QuadraticEigsolvOperators() {
00052 OutputCapturer capt;
00053 capt.begin_capture();
00054
00055 delete A;
00056 delete M;
00057 delete K;
00058 delete K11;
00059 delete K12;
00060 delete K22;
00061 delete Y;
00062 delete H;
00063 delete H11;
00064 delete H12;
00065 delete H22;
00066 delete C;
00067 delete h_solver_op_;
00068 delete HPrec;
00069 delete Prec;
00070 delete Aprec;
00071 delete AprecSolver;
00072
00073 capt.end_capture();
00074
00075 if (capt.get_stdout().size() > 0)
00076 rDebug("%s", capt.get_stdout().c_str());
00077 }
00078
00079 int QuadraticEigsolvOperators::buildPreconditioner(double sigma) {
00080
00081 assert(K == 0);
00082 K = build_Asigma(*A, *M, sigma);
00083
00084
00085 Aprec = new Epetra_LinearProblem();
00086 Aprec->SetOperator(K);
00087
00088 AprecSolver = new AztecOO(*Aprec);
00089
00090
00091 AprecSolver->SetAztecOption(AZ_precond, AZ_Neumann);
00092 AprecSolver->SetAztecOption(AZ_poly_ord, 10);
00093
00094
00095
00096
00097
00098 AprecSolver->SetAztecOption(AZ_output, AZ_none);
00099
00100
00101 AprecSolver->SetAztecOption(AZ_solver, AZ_gmres);
00102
00103
00104 int iterAztec = 2;
00105 Prec = new AztecOO_Operator(AprecSolver, iterAztec);
00106
00107 return(0);
00108 }
00109
00110 void QuadraticEigsolvOperators::matrixAssembly(FemaxMesh& mesh_mixed) {
00111
00112 NedelecMesh& mesh_nedelec = mesh_mixed.get_nedelec_mesh();
00113 LagrangeMesh& mesh_lagrange = mesh_mixed.get_lagrange_mesh();
00114
00115
00116 std::vector<int> owned_dofs;
00117 for (mesh::id_t ldof = 0; ldof < mesh_nedelec.get_num_my_dofs(0); ++ ldof) {
00118 mesh::id_t gdof = mesh_nedelec.gdof(ldof);
00119 if (mesh_nedelec.is_owned_gdof(gdof)) {
00120 mesh::id_t mgdof = mesh_nedelec.map_dof(ldof);
00121 if (mgdof != mesh::ID_NONE)
00122 owned_dofs.push_back(mgdof);
00123 }
00124 }
00125 Epetra_Map map_nedelec (-1, owned_dofs.size(), &owned_dofs[0], 0, Comm);
00126 rInfoAll("map_nedelec: owned_dofs.size()=%u", static_cast<unsigned int>(owned_dofs.size()));
00127
00128 owned_dofs.clear();
00129 for (mesh::id_t ldof = 0; ldof < mesh_lagrange.get_num_my_dofs(0); ++ ldof) {
00130 mesh::id_t gdof = mesh_lagrange.gdof(ldof);
00131 if (mesh_lagrange.is_owned_gdof(gdof)) {
00132 mesh::id_t mgdof = mesh_lagrange.map_dof(ldof);
00133 if (mgdof != mesh::ID_NONE)
00134 owned_dofs.push_back(mgdof);
00135 }
00136 }
00137 Epetra_Map map_lagrange (-1, owned_dofs.size(), &owned_dofs[0], 0, Comm);
00138
00139
00140
00141
00142 rInfo("Assemble finite element matrices...");
00143
00144
00145 log_mem_footprint(Comm, "before matrix assembly");
00146 pbe_start(60, "assembly of A and M");
00147 A = new Epetra_FECrsMatrix(Copy, map_nedelec, false);
00148 M = new Epetra_FECrsMatrix(Copy, map_nedelec, false);
00149
00150 mesh_nedelec.assembleAM(A, M);
00151
00152
00153
00154 A->FillComplete(map_nedelec, map_nedelec);
00155 M->FillComplete(map_nedelec, map_nedelec);
00156 pbe_stop(60);
00157 log_matrix_stats("Matrix A", A);
00158 log_matrix_stats("Matrix M", M);
00159 log_mem_footprint(Comm, "after assembly of A and M");
00160
00161
00162 pbe_start(61, "construction of Y");
00163 Y = new Epetra_FECrsMatrix(Copy, map_nedelec, false);
00164 mesh_mixed.constructY(*Y, map_lagrange, map_nedelec);
00165
00166 Epetra_FECrsMatrix* Y_transp = new Epetra_FECrsMatrix (Copy, map_lagrange, false);
00167 mesh_mixed.constructY_transp(*Y_transp, map_nedelec, map_lagrange);
00168
00169 pbe_stop(61);
00170 log_matrix_stats("Matrix Y", Y);
00171 log_matrix_stats("Matrix Y'", Y_transp);
00172 log_mem_footprint(Comm, "after constr of Y and Y'");
00173
00174
00175
00176
00177
00178 pbe_start(62, "construction of C and C^T");
00179 C = new Epetra_CrsMatrix (Copy, map_nedelec, 0);
00180 EpetraExt::MatrixMatrix::Multiply(*M, false, *Y, false, *C);
00181 pbe_stop(62);
00182 log_matrix_stats("Matrix C", C);
00183 log_mem_footprint(Comm, "after constr of C");
00184
00185
00186 pbe_start(63, "construction of H");
00187 H = new Epetra_CrsMatrix (Copy, map_lagrange, 0);
00188 #ifdef USE_PM_MATRIXMATRIX
00189 poor_man_matrix_matrix(*Y_transp, *C, *H);
00190 #else
00191 EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *C, false, *H);
00192 #endif
00193 pbe_stop(63);
00194 log_matrix_stats("Matrix H", H);
00195 log_mem_footprint(Comm, "after constr of H");
00196
00197 if (params_.get<bool>("export_mtx")) {
00198 pbe_start(130, "MatrixMarket export");
00199 export_Epetra_CrsMatrix(*A, 15, "A.mtx");
00200 export_Epetra_CrsMatrix(*M, 15, "M.mtx");
00201 export_Epetra_CrsMatrix(*Y, 2, "Y.mtx");
00202 export_Epetra_CrsMatrix(*Y_transp, 2, "Y_transp.mtx");
00203 export_Epetra_CrsMatrix(*C, 15, "C.mtx");
00204 export_Epetra_CrsMatrix(*H, 15, "H.mtx");
00205
00206 export_Epetra_Map(A->RowMap(), "A_rowmap.txt");
00207 export_Epetra_Map(H->RowMap(), "H_rowmap.txt");
00208
00209 pbe_stop(130);
00210 }
00211
00212 delete Y_transp;
00213 log_mem_footprint(Comm, "after deleting Y'");
00214 }
00215
00216 int QuadraticEigsolvOperators::buildHierarchicalPrec(double sigma) {
00217
00218
00219
00220 bool is_ml = true;
00221
00222 const Epetra_Map& map_nedelec = M->OperatorDomainMap();
00223 int sizeLocal = map_nedelec.NumMyElements();
00224
00225 int nBlockRows1 = 0;
00226 int nBlockRows2 = 0;
00227 int i = 0;
00228 int gid = 0;
00229
00230 for (i = 0; i < sizeLocal; ++i) {
00231 if (map_nedelec.GID(i) < noe)
00232 nBlockRows1++;
00233 else
00234 nBlockRows2++;
00235 }
00236
00237 int *myBlockRow1 = new int[nBlockRows1];
00238 int *myBlockRow2 = new int[nBlockRows2];
00239
00240 nBlockRows1 = 0;
00241 nBlockRows2 = 0;
00242 for (i = 0; i < sizeLocal; ++i) {
00243 gid = map_nedelec.GID(i);
00244 if (gid < noe)
00245 myBlockRow1[nBlockRows1++] = gid;
00246 else
00247 myBlockRow2[nBlockRows2++] = gid-noe;
00248 }
00249
00250 Epetra_Map map1(noe, nBlockRows1, myBlockRow1, 0, Comm);
00251 Epetra_Map map2(map_nedelec.NumGlobalElements() - noe,
00252 nBlockRows2, myBlockRow2, 0, Comm);
00253
00254 K11 = new Epetra_CrsMatrix(Copy, map1, 0);
00255 K12 = new Epetra_CrsMatrix(Copy, map1, 0);
00256 K22 = new Epetra_CrsMatrix(Copy, map2, 0);
00257 if (is_ml)
00258 Y11 = new Epetra_CrsMatrix(Copy, map1, 0);
00259
00260 int numEntries = 0;
00261 int *indices = 0;
00262 double *values = 0;
00263 int info = 0;
00264 int j = 0;
00265 int col = 0;
00266 double val = 0.0;
00267
00268 int numEntriesY = 0;
00269 int *indicesY = 0;
00270 double *valuesY = 0;
00271
00272 for (i = 0; i < sizeLocal; ++i) {
00273 gid = map_nedelec.GID(i);
00274
00275 if (gid < noe) {
00276
00277 assert(A->ExtractMyRowView(i, numEntries, values, indices) >= 0 );
00278 for (j = 0; j < numEntries; ++j) {
00279 col = A->GCID(indices[j]);
00280 if (col < noe) {
00281 assert(K11->InsertGlobalValues(gid, 1, values + j, &col) >=0);
00282 }
00283 else {
00284 col -= noe;
00285 assert(K12->InsertGlobalValues(gid, 1, values + j, &col) >=0);
00286 }
00287 }
00288
00289 assert(M->ExtractMyRowView(i, numEntries, values, indices) >= 0);
00290 for (j = 0; j < numEntries; ++j) {
00291 col = M->GCID(indices[j]);
00292 val = - sigma*values[j];
00293 if (col < noe) {
00294 assert((info = K11->SumIntoGlobalValues(gid, 1, &val, &col)) >= 0);
00295 if (info > 0)
00296 assert(K11->InsertGlobalValues(gid, 1, &val, &col) >= 0);
00297 }
00298 else {
00299 col -= noe;
00300 assert((info = K12->SumIntoGlobalValues(gid, 1, &val, &col)) >= 0);
00301 if (info > 0)
00302 assert(K12->InsertGlobalValues(gid, 1, &val, &col) >= 0);
00303 }
00304 }
00305
00306 if (is_ml) {
00307 assert(Y->ExtractMyRowView(i, numEntriesY, valuesY, indicesY) >= 0 );
00308 for (j = 0; j < numEntriesY; ++j) {
00309 col = Y->GCID(indicesY[j]);
00310 if (col < noe) {
00311 assert(Y11->InsertGlobalValues(gid, 1, valuesY + j, &col) >=0);
00312 }
00313 }
00314 }
00315
00316 }
00317
00318 else {
00319
00320 gid -= noe;
00321
00322 assert(A->ExtractMyRowView(i, numEntries, values, indices) >=0);
00323 for (j = 0; j < numEntries; ++j) {
00324 col = A->GCID(indices[j]);
00325
00326 if (col >= noe && A->MyGRID(col)) {
00327 col -= noe;
00328 assert(K22->InsertGlobalValues(gid, 1, values + j, &col) >=0);
00329 }
00330 }
00331
00332 assert(M->ExtractMyRowView(i, numEntries, values, indices) >= 0);
00333 for (j = 0; j < numEntries; ++j) {
00334 col = M->GCID(indices[j]);
00335
00336
00337 if (col >= noe && A->MyGRID(col)) {
00338 col -= noe;
00339 val = - sigma*values[j];
00340 assert((info = K22->SumIntoGlobalValues(gid, 1, &val, &col)) >= 0);
00341 if (info > 0)
00342 assert(K22->InsertGlobalValues(gid, 1, &val, &col) >= 0);
00343 }
00344 }
00345 }
00346
00347 }
00348
00349 K11->FillComplete(map1, map1);
00350 K12->FillComplete(map2, map1);
00351 K22->FillComplete(map2, map2);
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 if (is_ml) {
00404 if (H11 == 0)
00405 H11 = getH11();
00406 Y11->FillComplete(H11->RangeMap(), map1);
00407 }
00408 Prec = new HierarchicalBasisPrec(K11, Y11, H11, K12, K22, M->OperatorDomainMap(), M->OperatorRangeMap(), params_.sublist("asigma_precon"));
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 delete myBlockRow1;
00503 delete myBlockRow2;
00504 return(0);
00505 }
00506
00507 int QuadraticEigsolvOperators::buildHierarchicalHPrec() {
00508
00509
00510
00511
00512
00513
00514
00515 const Epetra_Map& map_lagrange = H->OperatorDomainMap();
00516 int sizeLocal = map_lagrange.NumMyElements();
00517
00518 int nBlockRows1 = 0;
00519 int nBlockRows2 = 0;
00520 int i = 0;
00521 int gid = 0;
00522
00523 for (i = 0; i < sizeLocal; ++i) {
00524 if (map_lagrange.GID(i) < non)
00525 nBlockRows1++;
00526 else
00527 nBlockRows2++;
00528 }
00529
00530 int *myBlockRow1 = new int[nBlockRows1];
00531 int *myBlockRow2 = new int[nBlockRows2];
00532
00533 nBlockRows1 = 0;
00534 nBlockRows2 = 0;
00535 for (i = 0; i < sizeLocal; ++i) {
00536 gid = map_lagrange.GID(i);
00537 if (gid < non)
00538 myBlockRow1[nBlockRows1++] = gid;
00539 else
00540 myBlockRow2[nBlockRows2++] = gid-non;
00541 }
00542
00543 Epetra_Map map1(non, nBlockRows1, myBlockRow1, 0, Comm);
00544 Epetra_Map map2(noe, nBlockRows2, myBlockRow2, 0, Comm);
00545
00546 bool create_h11 = (H11 == 0) ? true : false;
00547 if (create_h11)
00548 H11 = new Epetra_CrsMatrix(Copy, map1, 0);
00549 H12 = new Epetra_CrsMatrix(Copy, map1, 0);
00550 H22 = new Epetra_CrsMatrix(Copy, map2, 0);
00551
00552 int numEntries = 0;
00553 int *indices = 0;
00554 double *values = 0;
00555 int j = 0;
00556 int col = 0;
00557
00558 for (i = 0; i < sizeLocal; ++i) {
00559 gid = map_lagrange.GID(i);
00560 if (gid < non) {
00561 assert(H->ExtractMyRowView(i, numEntries, values, indices) >= 0 );
00562 for (j = 0; j < numEntries; ++j) {
00563 col = H->GCID(indices[j]);
00564 if (col < non) {
00565 if (create_h11)
00566 assert(H11->InsertGlobalValues(gid, 1, values + j, &col) >=0);
00567 }
00568 else {
00569 col -= non;
00570 assert(H12->InsertGlobalValues(gid, 1, values + j, &col) >=0);
00571 }
00572 }
00573 }
00574 else {
00575 gid -= non;
00576 assert(H->ExtractMyRowView(i, numEntries, values, indices) >=0);
00577 for (j = 0; j < numEntries; ++j) {
00578 col = H->GCID(indices[j]);
00579
00580
00581 if (col >= non && H->MyGRID(col)) {
00582 col -= non;
00583 assert(H22->InsertGlobalValues(gid, 1, values + j, &col) >=0);
00584 }
00585 }
00586 }
00587 }
00588
00589 if (create_h11)
00590 H11->FillComplete(map1, map1);
00591 H12->FillComplete(map2, map1);
00592 H22->FillComplete(map2, map2);
00593
00594 HPrec = new HierarchicalBasisPrec(H11, 0, 0, H12, H22, H->OperatorDomainMap(), H->OperatorRangeMap(), params_.sublist("h_precon"));
00595
00596 delete myBlockRow1;
00597 delete myBlockRow2;
00598 return(0);
00599 }
00600
00601
00602 Epetra_CrsMatrix* QuadraticEigsolvOperators::getH11() {
00603
00604 if (H11 == 0) {
00605 const Epetra_Map& map_lagrange = H->OperatorDomainMap();
00606 int sizeLocal = map_lagrange.NumMyElements();
00607
00608 int nBlockRows1 = 0;
00609 int i = 0;
00610 int gid = 0;
00611
00612 for (i = 0; i < sizeLocal; ++i) {
00613 if (map_lagrange.GID(i) < non)
00614 nBlockRows1++;
00615 }
00616
00617 int *myBlockRow1 = new int[nBlockRows1];
00618
00619 nBlockRows1 = 0;
00620 for (i = 0; i < sizeLocal; ++i) {
00621 gid = map_lagrange.GID(i);
00622 if (gid < non)
00623 myBlockRow1[nBlockRows1++] = gid;
00624 }
00625
00626 Epetra_Map map1(non, nBlockRows1, myBlockRow1, 0, Comm);
00627 H11 = new Epetra_CrsMatrix(Copy, map1, 0);
00628
00629
00630 int numEntries = 0;
00631 int *indices = 0;
00632 double *values = 0;
00633 int j = 0;
00634 int col = 0;
00635
00636 for (i = 0; i < sizeLocal; ++i) {
00637 gid = map_lagrange.GID(i);
00638 if (gid < non) {
00639 assert(H->ExtractMyRowView(i, numEntries, values, indices) >= 0 );
00640 for (j = 0; j < numEntries; ++j) {
00641 col = H->GCID(indices[j]);
00642 if (col < non) {
00643 assert(H11->InsertGlobalValues(gid, 1, values + j, &col) >=0);
00644 }
00645 }
00646 }
00647 }
00648
00649 H11->FillComplete(map1, map1);
00650 }
00651
00652 return(H11);
00653 }
00654
00655 const Epetra_Operator* QuadraticEigsolvOperators::getA() const {
00656 return(A);
00657 }
00658
00659 const Epetra_Operator* QuadraticEigsolvOperators::getM() const {
00660 return(M);
00661 }
00662
00663 const Epetra_Operator* QuadraticEigsolvOperators::getAsigmaPrec(double sigma) {
00664 if (Prec == 0) {
00665 Teuchos::ParameterList& sublist(params_.sublist("asigma_precon"));
00666 string asigma_precon_type = sublist.get<string>("type");
00667 rInfo("Calculate preconditioner for A-sigma*M...");
00668 pbe_start(70, "Preconditioner A-sigma*M construction");
00669 if (asigma_precon_type == "neumann") {
00670 buildPreconditioner(sigma);
00671 } else if (asigma_precon_type == "2level") {
00672 buildHierarchicalPrec(sigma);
00673 } else if (asigma_precon_type == "lu") {
00674
00675 assert(K == 0);
00676 K = build_Asigma(*A, *M, sigma);
00677 Prec = new DirectSolverOperator(*K, params_.sublist("asigma_precon").sublist("amesos_params"));
00678 } else if (asigma_precon_type == "diag") {
00679
00680
00681 Prec = new DiagonalPreconditioner(*K);
00682 } else {
00683 rExit("Unknown preconditioner for A - sigma*M: %s", asigma_precon_type.c_str());
00684 }
00685 pbe_stop(70);
00686 } else {
00687 assert(_sigma_K == sigma);
00688 }
00689 _sigma_K = sigma;
00690 return(Prec);
00691 }
00692
00693 const Epetra_Operator* QuadraticEigsolvOperators::getHSolver() {
00694 if (h_solver_op_ == 0) {
00695 rInfo("Calculate solver for H...");
00696 pbe_start(72, "Solver H construction");
00697 Teuchos::ParameterList& sublist(params_.sublist("h_solver"));
00698 string h_solver_type = sublist.get<string>("type");
00699 if (h_solver_type == "pcg") {
00700
00701
00702 int max_iter_cg = 100;
00703 double h_solver_tol = 1e-10;
00704 int verbosity = 0;
00705 BlockPCGSolver* pcg = new BlockPCGSolver(Comm, H, h_solver_tol, max_iter_cg, verbosity);
00706 pcg->setPreconditioner(const_cast<Epetra_Operator*>(getHPrec()));
00707 h_solver_op_ = pcg;
00708 } else if (h_solver_type == "lu") {
00709
00710 h_solver_op_ = new DirectSolverOperator(*H, params_.sublist("h_solver").sublist("amesos_params"));
00711 } else if (h_solver_type == "prec") {
00712
00713 h_solver_op_ = const_cast<Epetra_Operator*>(getHPrec());
00714 } else {
00715 rExit("Unknown solver for H: %s", h_solver_type.c_str());
00716 }
00717 pbe_stop(72);
00718 }
00719 return(h_solver_op_);
00720 }
00721
00722 const Epetra_Operator* QuadraticEigsolvOperators::getHPrec() {
00723 if (HPrec == 0) {
00724 rInfo("Calculate preconditioner for H...");
00725 pbe_start(71, "Preconditioner H construction");
00726 Teuchos::ParameterList& sublist(params_.sublist("h_precon"));
00727 string h_precon_type = sublist.get<string>("type");
00728 if (h_precon_type == "2level") {
00729
00730 buildHierarchicalHPrec();
00731 } else if (h_precon_type == "ml") {
00732
00733 OutputCapturer cap;
00734 cap.begin_capture();
00735 HPrec = new ML_Epetra::MultiLevelPreconditioner(*H, params_.sublist("h_precon").sublist("ml_params"), true);
00736 cap.end_capture();
00737 cap.log("ML output:", RLOG_CHANNEL("debug"));
00738 } else {
00739 rExit("Unknown preconditioner for H: %s", h_precon_type.c_str());
00740 }
00741 pbe_stop(71);
00742 }
00743 return(HPrec);
00744 }
00745
00746 Epetra_Operator* QuadraticEigsolvOperators::getH() {
00747 return(H);
00748 }
00749
00750 Epetra_Operator* QuadraticEigsolvOperators::getY() const {
00751 return(Y);
00752 }
00753
00754 Epetra_Operator* QuadraticEigsolvOperators::getC() const {
00755 return(C);
00756 }