examples/matmatmult_bug.cpp

Go to the documentation of this file.
00001 //
00002 // C++ Implementation: matmatmult_bug
00003 //
00004 // Description:
00005 //
00006 //
00007 // Author: Roman Geus <roman.geus@psi.ch>, (C) 2004
00008 //
00009 
00010 #include <iostream>
00011 #include <fstream>
00012 #include <cstdlib>
00013 #include <unistd.h>
00014 
00015 #include "Epetra_ConfigDefs.h"
00016 #ifdef HAVE_MPI
00017 #include <mpi.h>
00018 #include "Epetra_MpiComm.h"
00019 #else
00020 #include "Epetra_SerialComm.h"
00021 #endif
00022 #include "Epetra_CrsMatrix.h"
00023 #include "Epetra_Map.h"
00024 #include "EpetraExt_MatrixMatrix.h"
00025 
00026 #include "balancedmtxreader.h"
00027 #include "mtxreader.h"
00028 #include "export_Epetra_CrsMatrix.h"
00029 
00030 int main(int argc, char *argv[]) {
00031 #ifdef HAVE_MPI
00032     MPI_Init(&argc, &argv);
00033     Epetra_MpiComm comm(MPI_COMM_WORLD);
00034 #else
00035     Epetra_SerialComm comm;
00036 #endif
00037 
00038     if (comm.MyPID() == 0)
00039         cout << "running on " << comm.NumProc() << " processor(s)..." << endl;
00040 
00041     BaseMtxReader *reader = 0;
00042 
00043     // Y matrix
00044     string filename("Y.mtx");
00045     if (access(filename.c_str(), R_OK)) {
00046         cerr << "Unable to access Y matrix file.\n";
00047         exit(1);
00048     }
00049     reader = new BalancedMtxReader(filename, comm);
00050     Epetra_CrsMatrix* y = reader->read();
00051     delete reader;
00052     const Epetra_Map& short_map = y->DomainMap();
00053     const Epetra_Map& long_map = y->RangeMap();
00054 
00055     // Y_transp matrix
00056     filename = "Y_transp.mtx";
00057     if (access(filename.c_str(), R_OK)) {
00058         cerr << "Unable to access Y_transp matrix file.\n";
00059         exit(1);
00060     }
00061     reader = new MtxReader(filename, short_map, long_map, short_map);
00062     Epetra_CrsMatrix* y_transp = reader->read();
00063     delete reader;
00064 
00065     // C matrix (with compatible maps)
00066     filename = "C.mtx";
00067     if (access(filename.c_str(), R_OK)) {
00068         cerr << "Unable to access C matrix file.\n";
00069         exit(1);
00070     }
00071     reader = new MtxReader(filename, long_map, short_map, long_map);
00072     Epetra_CrsMatrix* c = reader->read();
00073     delete reader;
00074 
00075     Epetra_CrsMatrix* ha = new Epetra_CrsMatrix (Copy, short_map, 0);
00076     EpetraExt::MatrixMatrix::Multiply(*y, true, *c, false, *ha);
00077 
00078     Epetra_CrsMatrix* hb = new Epetra_CrsMatrix (Copy, short_map, 0);
00079     EpetraExt::MatrixMatrix::Multiply(*y_transp, false, *c, false, *hb);
00080 
00081     export_Epetra_CrsMatrix(*ha, 5, "ha.mtx");
00082     export_Epetra_CrsMatrix(*hb, 5, "hb.mtx");
00083 
00084     delete y;
00085     delete y_transp;
00086     delete c;
00087     delete ha;
00088     delete hb;
00089 
00090     return 0;
00091 }

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