FEMAXX (Finite Element Maxwell Eigensolver): test/test_matrix

00001 //
00002 // C++ Implementation: test_romberg
00003 //
00004 // Description: unit tests for matrix-matrix operations
00005 //
00006 //
00007 // Author: Roman Geus <geus@maxwell>, (C) 2004
00008 //
00009 // Copyright: See COPYING file that comes with this distribution
00010 //
00011 //
00012 
00013 #include <cmath>
00014 #include <iostream>
00015 #include <fstream>
00016 #include <cstdlib>
00017 #include <unistd.h>
00018 #include <limits>
00019 #include <tut.h>
00020 #include "rlog/rlog.h"
00021 #include "Epetra_Comm.h"
00022 #include "Epetra_CrsMatrix.h"
00023 #include "Epetra_Map.h"
00024 #include "Epetra_Vector.h"
00025 #include "EpetraExt_MatrixMatrix.h"
00026 #include "EpetraExt_Transpose_RowMatrix.h"
00027 #include "AbstractEigsolvOperators.h"
00028 #include "balancedmtxreader.h"
00029 #include "mtxreader.h"
00030 #include "matrix_matrix.h"
00031 
00032 #undef RLOG_TIME_TSC
00033 #include "rlog/RLogTime.h"
00034 
00035 namespace tut {
00036     Epetra_Comm& get_comm_world();
00037 }
00038 
00039 namespace
00040 {
00041 
00042     struct testdata
00043     {
00044         testdata()
00045             : comm(tut::get_comm_world()),
00046               tol(100 * std::numeric_limits<double>::epsilon())
00047         {
00048             if (M == 0)
00049                 read_matrices(comm);
00050         }
00051 
00052         bool is_same_matrix(const Epetra_Operator& A, const Epetra_Operator& B) {
00053             int info;
00054 
00055             // random vectors x and y
00056             Epetra_Vector y(A.OperatorDomainMap());
00057             Epetra_Vector x(A.OperatorRangeMap());
00058             Epetra_Vector t(A.OperatorRangeMap());
00059             x.Random();
00060             y.Random();
00061 
00062             double s1, s2;
00063             A.Apply(y,t);
00064             info = t.Dot(x, &s1);
00065             rAssert(info == 0);
00066         
00067             B.Apply(y,t);
00068             info = t.Dot(x, &s2);
00069             rAssert(info == 0);
00070             rDebug("is_same_matrix: s1=%e, s2=%e, |s1-s2|=%e", s1, s2, fabs(s1-s2));
00071             return (fabs(s1-s2) < tol);
00072         }
00073 
00074         bool is_same_matrix_transp(const Epetra_Operator& A, const Epetra_Operator& B) {
00075             int info;
00076 
00077             // random vectors x and y
00078             Epetra_Vector y(A.OperatorDomainMap());
00079             Epetra_Vector x(A.OperatorRangeMap());
00080             Epetra_Vector t1(A.OperatorRangeMap());
00081             Epetra_Vector t2(A.OperatorDomainMap());
00082             x.Random();
00083             y.Random();
00084 
00085             double s1, s2;
00086             A.Apply(y,t1);
00087             info = t1.Dot(x, &s1);
00088             rAssert(info == 0);
00089         
00090             B.Apply(x,t2);
00091             info = t2.Dot(y, &s2);
00092             rAssert(info == 0);
00093             rDebug("is_same_matrix_transp: s1=%e, s2=%e, |s1-s2|=%e", s1, s2, fabs(s1-s2));
00094             return (fabs(s1-s2) < tol);
00095         }
00096 
00097         static void read_matrices(Epetra_Comm& comm) {
00098             rInfo("Reading test matrices...");
00099             string filename;
00100             BaseMtxReader *reader = 0;
00101 
00102             // read matrix Y_transp
00103             filename = "test/matrix_matrix/Y_transp.mtx";
00104             if (access(filename.c_str(), R_OK)) {
00105                 rError("Unable to access Y_transp matrix file.");
00106                 exit(1);
00107             }
00108             reader = new BalancedMtxReader(filename, comm);
00109             Y_transp = reader->read();
00110             delete reader;
00111             const Epetra_Map& long_map = Y_transp->DomainMap();
00112             const Epetra_Map& short_map = Y_transp->RangeMap();
00113 
00114             // read matrix Y
00115             filename = "test/matrix_matrix/Y.mtx";
00116             if (access(filename.c_str(), R_OK)) {
00117                 rError("Unable to access Y matrix file.");
00118                 exit(1);
00119             }
00120             reader = new MtxReader(filename, long_map, short_map, long_map);
00121             Y = reader->read();
00122             delete reader;
00123 
00124             // read matrix M
00125             filename = "test/matrix_matrix/M.mtx";
00126             if (access(filename.c_str(), R_OK)) {
00127                 rError("Unable to access Y matrix file.");
00128                 exit(1);
00129             }
00130             reader = new MtxReader(filename, long_map, long_map, long_map);
00131             M = reader->read();
00132             delete reader;
00133         }
00134 
00135         Epetra_Comm& comm;
00136         const double tol;
00137         static Epetra_CrsMatrix* Y_transp;
00138         static Epetra_CrsMatrix* Y;
00139         static Epetra_CrsMatrix* M;
00140     };
00141 
00142     // Initialise static members
00143     Epetra_CrsMatrix* testdata::Y_transp = 0;
00144     Epetra_CrsMatrix* testdata::Y = 0;
00145     Epetra_CrsMatrix* testdata::M = 0;
00146 
00147     typedef tut::test_group<testdata> testgroup;
00148     typedef testgroup::object testobject;
00149     testgroup group("matrix_matrix");
00150 
00152     template<>
00153     template<>
00154     void testobject::test<11>()
00155     {
00156         const Epetra_Map& short_map = Y_transp->RangeMap();
00157         int info;
00158         rlog_time t1, t2;
00159 
00160 #if 0
00161         // compute C0
00162         rlog_get_time(&t1);
00163         Epetra_CrsMatrix C0(Copy, short_map, 0);
00164         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *M, true, C0);
00165         rAssert(info == 0);
00166         rlog_get_time(&t2);
00167         rDebug("MatrixMatrix::Multiply [F,T]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00168 #endif
00169         // compute C1
00170         rlog_get_time(&t1);
00171         Epetra_CrsMatrix C1(Copy, short_map, 0);
00172         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *M, false, C1);
00173         rAssert(info == 0);
00174         rlog_get_time(&t2);
00175         rDebug("MatrixMatrix::Multiply [F,F]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00176 
00177         // compute C2
00178         rlog_get_time(&t1);
00179         Epetra_CrsMatrix C2(Copy, short_map, 0);
00180         info = poor_man_matrix_matrix(*Y_transp, *M, C2);
00181         rAssert(info == 0);
00182         rlog_get_time(&t2);
00183         rDebug("poor_man_matrix_matrix: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00184 
00185         // Test C1 == C2
00186 #if 0
00187         ensure(is_same_matrix(C0, C1));
00188 #endif
00189         ensure(is_same_matrix(C1, C2));
00190     }
00191 
00193     template<>
00194     template<>
00195     void testobject::test<12>()
00196     {
00197         const Epetra_Map& long_map = Y_transp->DomainMap();
00198         const Epetra_Map& short_map = Y_transp->RangeMap();
00199         int info;
00200         rlog_time t1, t2;
00201 
00202         // compute C1
00203         rlog_get_time(&t1);
00204         Epetra_CrsMatrix C1(Copy, short_map, 0);
00205         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *M, false, C1);
00206         rAssert(info == 0);
00207         rlog_get_time(&t2);
00208 
00209         // compute C2
00210         rlog_get_time(&t1);
00211         Epetra_CrsMatrix C2(Copy, short_map, 0);
00212         info = poor_man_matrix_matrix_local(*Y_transp, *M, C2);
00213         rAssert(info == 0);
00214         rlog_get_time(&t2);
00215         rDebug("poor_man_matrix_matrix_local: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00216 
00217         // compute C3
00218         rlog_get_time(&t1);
00219         Epetra_CrsGraph G3(Copy, short_map, M->ColMap(), 0);
00220         info = poor_man_matrix_matrix_graph(Y_transp->Graph(), M->Graph(), G3);
00221         rAssert(info == 0);
00222         info = G3.FillComplete(long_map, short_map);
00223         rAssert(info == 0);
00224         Epetra_CrsMatrix C3(Copy, G3);
00225         info = poor_man_matrix_matrix_local2(*Y_transp, *M, C3);
00226         rAssert(info == 0);
00227         rlog_get_time(&t2);
00228         rDebug("poor_man_matrix_matrix_local2: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00229 
00230         // Test C1 == C3
00231         ensure(is_same_matrix(C2, C3));
00232         if (comm.NumProc() == 1)
00233             ensure(is_same_matrix(C1, C3));
00234     }
00235 
00237     template<>
00238     template<>
00239     void testobject::test<13>()
00240     {
00241         const Epetra_Map& long_map = Y_transp->DomainMap();
00242         const Epetra_Map& short_map = Y_transp->RangeMap();
00243         int info;
00244         rlog_time t1, t2;
00245 
00246         // compute C' = Y'*M
00247         rlog_get_time(&t1);
00248         Epetra_CrsMatrix Ct(Copy, short_map, 0);
00249         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *M, false, Ct);
00250         rAssert(info == 0);
00251         rlog_get_time(&t2);
00252         rDebug("MatrixMatrix::Multiply Y'*M: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00253         // compute C = M*Y
00254         rlog_get_time(&t1);
00255         Epetra_CrsMatrix C(Copy, long_map, 0);
00256         info = EpetraExt::MatrixMatrix::Multiply(*M, false, *Y, false, C);
00257         rAssert(info == 0);
00258         rlog_get_time(&t2);
00259         rDebug("MatrixMatrix::Multiply M*Y: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00260         // Test C' == (C)' using random vectors
00261         ensure(is_same_matrix_transp(C, Ct));
00262     }
00263 
00265     template<>
00266     template<>
00267     void testobject::test<14>()
00268     {
00269         const Epetra_Map& short_map = Y_transp->RangeMap();
00270         int info;
00271         rlog_time t1, t2;
00272 
00273         // compute C' = Y'*M non-transposed MMM
00274         rlog_get_time(&t1);
00275         Epetra_CrsMatrix Ct1(Copy, short_map, 0);
00276         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *M, false, Ct1);
00277         rAssert(info == 0);
00278         rlog_get_time(&t2);
00279         rDebug("MatrixMatrix::Multiply Y'*M [F,F]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00280         // compute C' = Y'*M transposed MMM
00281         rlog_get_time(&t1);
00282         Epetra_CrsMatrix Ct2(Copy, short_map, 0);
00283         info = EpetraExt::MatrixMatrix::Multiply(*Y, true, *M, false, Ct2);
00284         rAssert(info == 0);
00285         rlog_get_time(&t2);
00286         rDebug("MatrixMatrix::Multiply Y'*M [T,F]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00287         // Test Ct1 == Ct2 using random vectors
00288         ensure(is_same_matrix(Ct1, Ct2));
00289     }
00290 
00292     template<>
00293     template<>
00294     void testobject::test<20>()
00295     {
00296         const Epetra_Map& long_map = Y_transp->DomainMap();
00297         const Epetra_Map& short_map = Y_transp->RangeMap();
00298         int info;
00299         rlog_time t1, t2;
00300 
00301         // compute C = M*Y
00302         Epetra_CrsMatrix C(Copy, long_map, 0);
00303         info = EpetraExt::MatrixMatrix::Multiply(*M, false, *Y, false, C);
00304         rAssert(info == 0);
00305         // compute Ct = Yt*M
00306         Epetra_CrsMatrix Ct(Copy, short_map, 0);
00307         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *M, false, Ct);
00308         rAssert(info == 0);
00309         // compute H1 = Y'*C non-transposed MMM
00310         rlog_get_time(&t1);
00311         Epetra_CrsMatrix H1(Copy, short_map, 0);
00312         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, C, false, H1);
00313         rAssert(info == 0);
00314         rlog_get_time(&t2);
00315         rDebug("MatrixMatrix::Multiply Y'*C [F,F]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00316         // compute H2 = Y'*C transposed MMM
00317         rlog_get_time(&t1);
00318         Epetra_CrsMatrix H2(Copy, short_map, 0);
00319         info = EpetraExt::MatrixMatrix::Multiply(*Y, true, C, false, H2);
00320         rAssert(info == 0);
00321         rlog_get_time(&t2);
00322         rDebug("MatrixMatrix::Multiply Y'*C [T,F]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00323         // compute H3 = C'*Y non-transposed MMM
00324         rlog_get_time(&t1);
00325         Epetra_CrsMatrix H3(Copy, short_map, 0);
00326         info = EpetraExt::MatrixMatrix::Multiply(Ct, false, *Y, false, H3);
00327         rAssert(info == 0);
00328         rlog_get_time(&t2);
00329         rDebug("MatrixMatrix::Multiply C'*Y [F,F]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00330         // compute H4 = C'*Y transposed MMM
00331         rlog_get_time(&t1);
00332         Epetra_CrsMatrix H4(Copy, short_map, 0);
00333         info = EpetraExt::MatrixMatrix::Multiply(C, true, *Y, false, H4);
00334         rAssert(info == 0);
00335         rlog_get_time(&t2);
00336         rDebug("MatrixMatrix::Multiply C'*Y [T,F]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00337         // Test H1 == H2 == H3 == H4 using random vectors
00338         ensure(is_same_matrix(H1, H2));
00339         ensure(is_same_matrix(H1, H3));
00340         ensure(is_same_matrix(H1, H4));
00341     }
00342 
00343 
00345     template<>
00346     template<>
00347     void testobject::test<22>()
00348     {
00349         const Epetra_Map& long_map = Y_transp->DomainMap();
00350         const Epetra_Map& short_map = Y_transp->RangeMap();
00351         int info;
00352         rlog_time t1, t2;
00353 
00354         // Compute C
00355         Epetra_CrsMatrix C(Copy, long_map, 0);
00356         info = EpetraExt::MatrixMatrix::Multiply(*M, false, *Y, false, C);
00357         rAssert(info == 0);
00358 
00359         // compute H1
00360         rlog_get_time(&t1);
00361         Epetra_CrsMatrix H1(Copy, short_map, 0);
00362         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, C, false, H1);
00363         rAssert(info == 0);
00364         rlog_get_time(&t2);
00365         rDebug("MatrixMatrix::Multiply [F,F]: H=Y'*C time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00366 
00367         // compute H2
00368         rlog_get_time(&t1);
00369         Epetra_CrsMatrix H2(Copy, short_map, 0);
00370         info = poor_man_matrix_matrix(*Y_transp, C, H2);
00371         rAssert(info == 0);
00372         rlog_get_time(&t2);
00373         rDebug("poor_man_matrix_matrix: H=Y'*C time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00374 
00375         // Test H1 == H2
00376         ensure(is_same_matrix(H1, H2));
00377     }
00378 
00380     template<>
00381     template<>
00382     void testobject::test<30>()
00383     {
00384         int info;
00385         rlog_time t1, t2;
00386 
00387         const bool MakeDataContiguous = true;
00388         const bool IgnoreNonLocalCols = false;
00389 
00390         const Epetra_Map& long_map = Y_transp->DomainMap();
00391         const Epetra_Map& short_map = Y_transp->RangeMap();
00392 
00393         Epetra_CrsMatrix Ct(Copy, short_map, 0);
00394         info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *M, false, Ct);
00395         rAssert(info == 0);
00396 
00397         Epetra_CrsMatrix C(Copy, long_map, 0);
00398         info = EpetraExt::MatrixMatrix::Multiply(*M, false, *Y, false, C);
00399         rAssert(info == 0);
00400 
00401         rlog_get_time(&t1);
00402         Epetra_Map map(Y->RowMap());
00403         EpetraExt::RowMatrix_Transpose transposer(MakeDataContiguous,
00404                                                   &map,
00405                                                   IgnoreNonLocalCols);
00406         Epetra_RowMatrix& C_(transposer(Ct));
00407         rlog_get_time(&t2);
00408         rDebug("RowMatrix_Transpose C' [T,F]: time=%ld " RLOG_TIME_UNIT, rlog_time_diff(t2, t1));
00409         ensure(is_same_matrix(C, C_));
00410         ensure(is_same_matrix_transp(Ct, C_));
00411     }
00412 
00413 };
test/test_matrix_matrix.cpp
