src/eigsolv/AbstractEigsolvOperators.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Implementation: AbstractEigsolvOperators
00003 //
00004 // Description:
00005 //
00006 //
00007 // Author: Roman Geus <roman.geus@psi.ch>, (C) 2004
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
00010 //
00011 //
00012 
00013 #include <cassert>
00014 #include "myrlog.h"
00015 #include "Epetra_Comm.h"
00016 #include "Epetra_CrsMatrix.h"
00017 #include "Epetra_Map.h"
00018 #include "AbstractEigsolvOperators.h"
00019 
00020 #if 1
00021 
00022 Epetra_CrsMatrix* AbstractEigsolvOperators::build_Asigma(const Epetra_CrsMatrix& A,
00023                                                          const Epetra_CrsMatrix& M,
00024                                                          double sigma)
00025 {
00026     // FIXME: This routine could probably be sped up by inserting/adding whole rows instead of
00027     //        individual elements.
00028     const Epetra_Map& map_nedelec = A.DomainMap();
00029     Epetra_CrsMatrix* Asigma = new Epetra_CrsMatrix(Copy, map_nedelec, 0);
00030     int sizeLocal = map_nedelec.NumMyElements();
00031     int i;
00032     for (i = 0; i < sizeLocal; ++ i) {
00033         int numEntries = 0;
00034         int *indices;
00035         double *values;
00036         int i_global = map_nedelec.GID(i);
00037         int info;
00038         info = A.ExtractMyRowView(i, numEntries, values, indices);
00039         assert(info >= 0);
00040         int jj;
00041         for (jj = 0; jj < numEntries; ++ jj) {
00042             int col = A.GCID(indices[jj]);
00043             assert(Asigma->InsertGlobalValues(i_global, 1, values + jj, &col) >=0);
00044         }
00045         info = M.ExtractMyRowView(i, numEntries, values, indices);
00046         assert(info >= 0);
00047         for (jj = 0; jj < numEntries; ++ jj) {
00048             int col = M.GCID(indices[jj]);
00049             double val = - sigma*values[jj];
00050             info = Asigma->SumIntoGlobalValues(i_global, 1, &val, &col);
00051             assert(info >= 0);
00052             if (info > 0) {
00053                 assert(Asigma->InsertGlobalValues(i_global, 1, &val, &col) >= 0);
00054             }
00055         }
00056     }
00057 
00058     Asigma->FillComplete(map_nedelec, map_nedelec);
00059     log_matrix_stats("Matrix K", Asigma);
00060     return Asigma;
00061 }
00062 #else
00063 
00064 #include "EpetraExt_MatrixMatrix.h"
00065 
00066 Epetra_CrsMatrix* AbstractEigsolvOperators::build_Asigma(const Epetra_CrsMatrix& A,
00067                                                          const Epetra_CrsMatrix& M,
00068                                                          double sigma)
00069 {
00070     Epetra_CrsMatrix* Asigma = new Epetra_CrsMatrix(M);
00071     assert(EpetraExt::MatrixMatrix::Add(A, false, 1.0, *Asigma, -sigma) == 0);
00072     return Asigma;
00073 }
00074 
00075 #endif

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