src/eigsolv/QuadraticEigsolvOperators.cpp

Go to the documentation of this file.
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     // Logging empty strings causes a "[FMT ERROR]: failed to format msg" message in rlog
00075     if (capt.get_stdout().size() > 0)
00076         rDebug("%s", capt.get_stdout().c_str());
00077 }
00078 
00079 int QuadraticEigsolvOperators::buildPreconditioner(double sigma) {
00080     // construct K = A - sigma*M
00081     assert(K == 0);
00082     K = build_Asigma(*A, *M, sigma);
00083 
00084     // Setting the AZ_Neumann preconditioner
00085     Aprec = new Epetra_LinearProblem();
00086     Aprec->SetOperator(K);
00087 
00088     AprecSolver = new AztecOO(*Aprec);
00089 
00090     // Specify the preconditioner inside AztecOO
00091     AprecSolver->SetAztecOption(AZ_precond, AZ_Neumann);
00092     AprecSolver->SetAztecOption(AZ_poly_ord, 10);
00093 
00094     //  AprecSolver.SetAztecOption(AZ_precond, AZ_dom_decomp);
00095     //  AprecSolver.SetAztecOption(AZ_ilu, AZ_subdomain_solve);
00096 
00097     // Specify the output level
00098     AprecSolver->SetAztecOption(AZ_output, AZ_none);
00099 
00100     // Specify the iterative solver used
00101     AprecSolver->SetAztecOption(AZ_solver, AZ_gmres);
00102 
00103     // Define an operator for this preconditioner
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     // Create Epetra_Maps
00116     std::vector<int> owned_dofs;    // Epetra_Map expects int array
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     // Matrix assembly
00141     //
00142     rInfo("Assemble finite element matrices...");
00143 
00144     // Allocate matrices A and M
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     // Finish data entry
00153     // transforms to local index space, sorts column indices ans merges redundant entries
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     // Calculate Y
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     // Passing true to arg2 of MatrixMatrix::Multiply should be avoided since it slows down the routine by a factor of ~10
00175     // This is the reason why Y_transp is constructed.
00176 
00177     // Calculate C = M * Y
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     // Calculate H = Y^T * C
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     // We conservatively assume, that ML will be used to solve the (1,1)-block. In 
00218     // doing so we guarantee, that Y11 and H11 are constructed. This will make sure that 
00219     // Maxwell ML will be used in HierarchicalBasisPrec.
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                 //              if (col >= noe) {
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                 //              if (col >= noe) {
00336                 // FIXME: Only local diagonal blocks stored in K22!
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       int numProc = Comm.NumProc();
00355       int myPid = Comm.MyPID();
00356       int* sizes = new int[numProc];
00357       sizes[myPid] = map2.NumMyElements();
00358       Comm.GatherAll(sizes+myPid, sizes, 1);
00359       int* delim = new int[numProc+1];
00360       delim[0] = 0;
00361       for (i=0; i<numProc; i++)
00362       delim[i+1] = delim[i]+sizes[i];
00363 
00364       if ( myPid==0 ) {
00365       cout << "Delimiters : ";
00366       for (i=0; i<=numProc; i++)
00367       cout << delim[i] << " ";
00368       cout << endl;
00369       }
00370       Comm.Barrier();
00371 
00372       int k = 0;
00373       int diag = 0;
00374       int numDiag = 0;
00375       int numOutDiag = 0;
00376       int in = 0;
00377       int out = 0;
00378       for (i=0; i<map2.NumMyElements(); i++) {
00379       gid = map2.GID(i);
00380 
00381       for (diag=1; gid >= delim[diag]; diag++);
00382       K22->ExtractMyRowView(i, numEntries, values, indices);
00383       for (j = 0; j < numEntries; ++j) {
00384       col = K22->GCID(indices[j]);
00385       for (k=1; col >= delim[k]; k++);
00386       if (diag == k) numDiag++; else numOutDiag++;
00387       }
00388 
00389       if ( gid >= delim[myPid] && gid < delim[myPid+1] )  
00390       in++;
00391       else
00392       out++;
00393       }
00394 
00395       cout << myPid << "  in = " << in << "   out = " << out << endl;
00396 
00397       cout << myPid << " numDiag = " << numDiag << "  numOutDiag = " << numOutDiag << endl;
00398 
00399       delete delim;
00400       delete sizes;
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       K = new Epetra_CrsMatrix(Copy, map_nedelec, 0);
00413       sizeLocal = map_nedelec.NumMyElements();
00414       for (i = 0; i < sizeLocal; ++i) {
00415       int numEntries = 0;
00416       int *indices;
00417       double *values;
00418       int iMap = map_nedelec.GID(i);
00419       int info;
00420       info = A->ExtractMyRowView(i, numEntries, values, indices);
00421       assert(info >= 0);
00422       int jj;
00423       for (jj = 0; jj < numEntries; ++jj) {
00424       int col = A->GCID(indices[jj]);
00425       assert(K->InsertGlobalValues(iMap, 1, values + jj, &col) >=0);
00426       }
00427       info = M->ExtractMyRowView(i, numEntries, values, indices);
00428       assert(info >= 0);
00429       for (jj = 0; jj < numEntries; ++jj) {
00430       int col = M->GCID(indices[jj]);
00431       double val = - sigma*values[jj];
00432       info = K->SumIntoGlobalValues(iMap, 1, &val, &col);
00433       assert(info >= 0);
00434       if (info > 0) {
00435       assert(K->InsertGlobalValues(iMap, 1, &val, &col) >= 0);
00436       }
00437       }
00438       }
00439       K->FillComplete(map_nedelec, map_nedelec);
00440     */
00441 
00442     /*
00443 
00444     if (Comm.MyPID() == 0) {
00445     cout << "TEST K :......................" << endl;
00446     }
00447     Epetra_MultiVector X(map_nedelec, 1, false);
00448     Epetra_MultiVector Y1(map_nedelec, 1, false);
00449     Epetra_MultiVector Y2(map_nedelec, 1, false);
00450     X.Random();
00451 
00452 
00453     Epetra_LinearProblem Problem1;
00454     Problem1.SetOperator( K );
00455     Problem1.SetLHS( &Y1 );
00456     Problem1.SetRHS( &X );
00457 
00458     AztecOO Solver(Problem1);
00459     Solver.SetAztecOption( AZ_solver,  AZ_gmres );
00460     //    Solver.SetAztecOption( AZ_output, AZ_none );
00461     Solver.Iterate(100, 1e-10);
00462 
00463     if (Comm.MyPID() == 0) {
00464     cout << "TEST Prec :......................" << endl;
00465     }
00466 
00467     Prec->ApplyInverse(X, Y2);
00468 
00469     if (Comm.MyPID() == 0) {
00470     cout << "TEST Prec done :......................" << endl;
00471     }
00472 
00473     double nrm1 = 0.0;
00474     double nrm2 = 0.0;
00475     Y1.Norm2(&nrm1);
00476     Y2.Norm2(&nrm2);
00477     if (Comm.MyPID() == 0) {
00478     cout << "Norm(Y1) = " << nrm1 << endl;
00479     cout << "Norm(Y2) = " << nrm2 << endl;
00480     }    
00481     exit(0);
00482     */
00483 
00484     /*
00485       Epetra_Vector X(map_nedelec, false);
00486       Epetra_Vector Y1(map_nedelec, false);
00487       Epetra_Vector Y2(map_nedelec, false);
00488       X.Random();
00489       K->Multiply(false, X, Y1);
00490       Prec->Apply(X, Y2);
00491       double nrm1 = 0.0;
00492       double nrm2 = 0.0;
00493       Y1.Norm2(&nrm1);
00494       Y2.Norm2(&nrm2);
00495       if (Comm.MyPID() == 0) {
00496       cout << "Norm(Y1) = " << nrm1 << endl;
00497       cout << "Norm(Y2) = " << nrm2 << endl;
00498       }    
00499       exit(0);
00500     */    
00501 
00502     delete myBlockRow1;
00503     delete myBlockRow2;
00504     return(0);
00505 }
00506 
00507 int QuadraticEigsolvOperators::buildHierarchicalHPrec() {
00508 
00509     //
00510     // construct H = | H11 H12 |
00511     //               | H12 H22 |
00512     //
00513     // dim(H11) = non    dim(H22) = noe 
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                 //              if (col >= non) {
00580                 // FIXME: Only local diagonal blocks stored in H22!
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             // Preconditioner is LU-factorisation of K = A -sigma*M
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             // Preconditioner is diagonal of K = A -sigma*M
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             // block PCG solver
00701             // Fixme: parametrize this....
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             // direct solver
00710             h_solver_op_ = new DirectSolverOperator(*H, params_.sublist("h_solver").sublist("amesos_params"));
00711         } else if (h_solver_type == "prec") {
00712             // just do one preconditioning step instead of a full solve
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             // 2-level hierarchical preconditioner
00730             buildHierarchicalHPrec();
00731         } else if (h_precon_type == "ml") {
00732             // ML preconditioner for positive-definite matrix H
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 }

Generated on Fri Oct 26 13:35:11 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7