00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
00302 Epetra_CrsMatrix C(Copy, long_map, 0);
00303 info = EpetraExt::MatrixMatrix::Multiply(*M, false, *Y, false, C);
00304 rAssert(info == 0);
00305
00306 Epetra_CrsMatrix Ct(Copy, short_map, 0);
00307 info = EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *M, false, Ct);
00308 rAssert(info == 0);
00309
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
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
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
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
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
00355 Epetra_CrsMatrix C(Copy, long_map, 0);
00356 info = EpetraExt::MatrixMatrix::Multiply(*M, false, *Y, false, C);
00357 rAssert(info == 0);
00358
00359
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
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
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 };