00001
00002
00003
00004
00005
00006
00007
00008
00009
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
00027
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