EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/MatrixMatrix_LL/cxx_main.cpp
Go to the documentation of this file.
1/*
2//@HEADER
3// ***********************************************************************
4//
5// EpetraExt: Epetra Extended - Linear Algebra Services Package
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ***********************************************************************
41//@HEADER
42*/
43
44
45#include <string>
46#include <vector>
47#include <iostream>
48#include <fstream>
49
50#include <Epetra_ConfigDefs.h>
51
52#ifdef EPETRA_MPI
53#include <mpi.h>
54#include <Epetra_MpiComm.h>
55#endif
56
57#include <Epetra_SerialComm.h>
58#include <Epetra_Time.h>
59#include <Epetra_Import.h>
60#include <Epetra_Export.h>
61#include <Epetra_Map.h>
62#include <Epetra_LocalMap.h>
63#include <Epetra_CrsGraph.h>
64#include <Epetra_CrsMatrix.h>
65#include <Epetra_Vector.h>
66#include <EpetraExt_MMHelpers.h>
68
72
73namespace EpetraExt {
74extern
76 const Epetra_Map& column_map);
77}
78
80 const std::string& input_file_name,
81 std::vector<std::string>& filenames);
82
84 const std::string& input_file_name,
85 std::string& A_file,
86 bool& transA,
87 std::string& B_file,
88 bool& transB,
89 std::string& C_file);
90
91int broadcast_name(Epetra_Comm& Comm, std::string& name);
92
93int create_maps(Epetra_Comm& Comm,
94 const std::string& input_file_name,
95 Epetra_Map*& row_map,
96 Epetra_Map*& col_map,
97 Epetra_Map*& range_map,
98 Epetra_Map*& domain_map);
99
100int read_matrix(const std::string& filename,
101 Epetra_Comm& Comm,
102 const Epetra_Map* rowmap,
103 Epetra_Map* colmap,
104 const Epetra_Map* rangemap,
105 const Epetra_Map* domainmap,
106 Epetra_CrsMatrix*& mat);
107
108int run_test(Epetra_Comm& Comm,
109 const std::string& filename,
110 bool do_FillComplete,
111 bool result_mtx_to_file=false,
112 bool verbose=false);
113
114int two_proc_test(Epetra_Comm& Comm,
115 bool verbose=false);
116
117int test_find_rows(Epetra_Comm& Comm);
118
120 int localProc,
121 int local_n,
122 bool callFillComplete = true,
123 bool symmetric = true);
124
126 bool verbose);
127
128int test_drumm1(Epetra_Comm& Comm);
129
131//Global variable!!!!
132std::string path;
134
135int main(int argc, char** argv) {
136
137#ifdef EPETRA_MPI
138 MPI_Init(&argc,&argv);
139 Epetra_MpiComm Comm(MPI_COMM_WORLD);
140#else
142#endif
143
144 bool write_result_mtx = false;
145 bool verbose = false;
146 int write = 0, inp_specified = 0;
147 bool path_specified = false;
148 std::string input_file;
149 bool input_file_specified = false;
150
151 if (Comm.MyPID()==0) {
152 for(int ii=0; ii<argc; ++ii) {
153 if (!strcmp("-write_result", argv[ii])) write_result_mtx = true;
154 if (!strcmp("-v", argv[ii])) verbose = true;
155 if (!strcmp("-i", argv[ii])) {
156 input_file = argv[ii+1];
157 input_file_specified = true;
158 }
159 if (!strcmp("-d",argv[ii])) {
160 path = argv[ii+1];
161 path_specified = true;
162 }
163 }
164 write = write_result_mtx ? 1 : 0;
165 inp_specified = input_file_specified ? 1 : 0;
166 }
167#ifdef EPETRA_MPI
168 MPI_Bcast(&write, 1, MPI_INT, 0, MPI_COMM_WORLD);
169 if (write) write_result_mtx = true;
170 MPI_Bcast(&inp_specified, 1, MPI_INT, 0, MPI_COMM_WORLD);
171 if (inp_specified) input_file_specified = true;
172#endif
173
174 if (!path_specified) {
175 path = ".";
176 }
177
178 int err = two_proc_test(Comm, verbose);
179 if (err != 0) {
180 std::cout << "two_proc_test returned err=="<<err<<std::endl;
181 return(err);
182 }
183
184 std::vector<std::string> filenames;
185
186 if (input_file_specified) {
187 filenames.push_back(std::string(input_file));
188 }
189 else {
190 input_file = "infiles";
191 std::cout << "specific input-file not specified, getting list of input-files from file 'infiles'." << std::endl;
192
193 err = read_input_file(Comm, input_file, filenames);
194 if (err != 0) {
195 if (path_specified) path_specified = false;
196 path = "./MatrixMatrix";
197 read_input_file(Comm, input_file, filenames);
198 }
199 }
200
201 err = test_find_rows(Comm);
202 if (err != 0) {
203 std::cout << "test_find_rows returned err=="<<err<<std::endl;
204 return(err);
205 }
206
207 err = test_drumm1(Comm);
208 if (err != 0) {
209 std::cout << "test_drumm1 returned err=="<<err<<std::endl;
210 return(err);
211 }
212
213 for(size_t i=0; i<filenames.size(); ++i) {
214 err = run_test(Comm, filenames[i], true, write_result_mtx, verbose);
215 if (err != 0) break;
216 err = run_test(Comm, filenames[i], false, write_result_mtx, verbose);
217 if (err != 0) break;
218 }
219
221 Comm.MyPID(), 10,
222 true, false);
223
224// std::cout << "D: \n" << *D << std::endl;
225
226 EpetraExt::MatrixMatrix::Add(*D, true, 0.5, *D, 0.5);
227
228// std::cout << "symm D: \n" << *D << std::endl;
229
230 delete D;
231
232 if (err == 0) {
233 err = time_matrix_matrix_multiply(Comm, verbose);
234 }
235
236 int global_err = err;
237
238#ifdef EPETRA_MPI
239 MPI_Allreduce(&err, &global_err, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
240 MPI_Finalize();
241#endif
242
243 return(global_err);
244}
245
247{
248 int numprocs = Comm.NumProc();
249 int localproc = Comm.MyPID();
250 int numlocalrows = 2;
251 long long numglobalrows = ((long long)numprocs)*numlocalrows;
252 Epetra_Map rowmap(numglobalrows, 0LL, Comm);
253 Epetra_CrsMatrix matrix(Copy, rowmap, numglobalrows);
254
255 int err = 0;
256 long long* cols = new long long[numglobalrows];
257 double*vals = new double[numglobalrows];
258
259 for(long long j=0; j<numglobalrows; ++j) {
260 cols[j] = j;
261 vals[j] = 1.0;
262 }
263
264 Epetra_Map colmap(-1LL, numglobalrows, cols, 0LL, Comm);
265
266 for(int i=0; i<numlocalrows; ++i) {
267 long long row = ((long long)localproc)*numlocalrows+i;
268 err = matrix.InsertGlobalValues(row, 1, &(vals[i]), &row);
269 if (err != 0) {
270 return(err);
271 }
272 }
273
274 err = matrix.FillComplete();
275 if (err != 0) {
276 return(err);
277 }
278
279 Epetra_Map* map_rows = EpetraExt::find_rows_containing_cols(matrix, colmap);
280
281 if (map_rows->NumMyElements() != numglobalrows) {
282 return(-1);
283 }
284
285 delete map_rows;
286 delete [] cols;
287 delete [] vals;
288
289 return(0);
290}
291
292int expand_name_list(const char* newname,
293 const char**& names,
294 int& alloc_len,
295 int& num_names)
296{
297 int offset = num_names;
298 if (offset >= alloc_len) {
299 int alloc_increment = 8;
300 const char** newlist = new const char*[alloc_len+alloc_increment];
301 for(int i=0; i<offset; ++i) {
302 newlist[i] = names[i];
303 }
304 delete [] names;
305 names = newlist;
306 alloc_len += alloc_increment;
307 for(int i=offset; i<alloc_len; ++i) {
308 names[i] = NULL;
309 }
310 }
311
312 names[offset] = newname;
313 ++num_names;
314 return(0);
315}
316
317int broadcast_name(Epetra_Comm& Comm, std::string& name)
318{
319 if (Comm.NumProc() < 2) return(0);
320
321#ifdef EPETRA_MPI
322 int len = name.size();
323 int localProc = Comm.MyPID();
324 if (localProc == 0) {
325 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
326 MPI_Bcast((void*)name.c_str(), len+1, MPI_CHAR, 0, MPI_COMM_WORLD);
327 }
328 else {
329 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
330 name.resize(len+1, ' ');
331 MPI_Bcast((void*)name.c_str(), len+1, MPI_CHAR, 0, MPI_COMM_WORLD);
332 }
333
334#endif
335 return(0);
336}
337
339 const std::string& input_file_name,
340 std::vector<std::string>& filenames)
341{
342 int local_err = 0, global_err = 0;
343 std::ifstream* infile = NULL;
344 int pathlen = path.size();
345
346 if (Comm.MyPID() == 0) {
347 std::string full_name = path.empty() ? input_file_name : path+"/"+input_file_name;
348
349 infile = new std::ifstream(full_name.c_str());
350 if (!(*infile)) {
351 local_err = -1;
352 delete infile;
353 }
354 }
355
356 Comm.SumAll(&local_err, &global_err, 1);
357
358 if (global_err != 0) {
359 return(global_err);
360 }
361
362
363 if (Comm.MyPID() == 0) {
364 const int linelen = 512;
365 char line[linelen];
366
367 std::ifstream& ifile = *infile;
368 while(!ifile.eof()) {
369 if (pathlen>0) {
370 sprintf(line,"%s/",path.c_str());
371 ifile.getline(&(line[pathlen+1]), linelen);
372 }
373 else {
374 ifile.getline(line, linelen);
375 }
376
377 if (ifile.fail()) {
378 break;
379 }
380 if (strchr(line, '#') == NULL) {
381 filenames.push_back(std::string(line));
382 }
383 }
384
385 int numfiles = filenames.size();
386#ifdef EPETRA_MPI
387 MPI_Bcast(&numfiles, 1, MPI_INT, 0, MPI_COMM_WORLD);
388#endif
389 for(int i=0; i<numfiles; ++i) {
390 broadcast_name(Comm, filenames[i]);
391 }
392
393 delete infile;
394 }
395 else {
396 int numfiles = 0;
397#ifdef EPETRA_MPI
398 MPI_Bcast(&numfiles, 1, MPI_INT, 0, MPI_COMM_WORLD);
399#endif
400 filenames.resize(numfiles);
401 for(int i=0; i<numfiles; ++i) {
402 broadcast_name(Comm, filenames[i]);
403 }
404 }
405
406 return(0);
407}
408
410 const std::string& filename,
411 bool do_FillComplete,
412 bool result_mtx_to_file,
413 bool verbose)
414{
415 std::string A_file;
416 char AT[3]; AT[0] = '^'; AT[1] = 'T'; AT[2] = '\0';
417 std::string B_file;
418 char BT[3]; BT[0] = '^'; BT[1] = 'T'; BT[2] = '\0';
419 std::string C_file;
420 bool transA, transB;
421
422 if(!Comm.MyPID()) std::cout<<"Testing: "<<filename<<std::endl;
423
424 int err = read_matrix_file_names(Comm, filename, A_file, transA,
425 B_file, transB, C_file);
426 if (err != 0) {
427 std::cout << "Error, read_matrix_file_names returned " << err << std::endl;
428 return(err);
429 }
430
431 if (!transA) AT[0] = '\0';
432 if (!transB) BT[0] = '\0';
433
434 int localProc = Comm.MyPID();
435
436 if (localProc == 0 && verbose) {
437 std::cout << "Testing C=A"<<AT<<"*B"<<BT<< "; A:" << A_file
438 << ", B:" << B_file << ", C:" << C_file << std::endl;
439 }
440
441 Epetra_CrsMatrix* A = NULL;
442 Epetra_CrsMatrix* B = NULL;
443 Epetra_CrsMatrix* C = NULL;
444 Epetra_CrsMatrix* C_check = NULL;
445
446 Epetra_Map* A_row_map = NULL;
447 Epetra_Map* A_col_map = NULL;
448 Epetra_Map* A_range_map = NULL;
449 Epetra_Map* A_domain_map = NULL;
450 err = create_maps(Comm, A_file, A_row_map, A_col_map, A_range_map, A_domain_map);
451 if (err != 0) {
452 std::cout << "create_maps A returned err=="<<err<<std::endl;
453 return(err);
454 }
455
456 Epetra_Map* B_row_map = NULL;
457 Epetra_Map* B_col_map = NULL;
458 Epetra_Map* B_range_map = NULL;
459 Epetra_Map* B_domain_map = NULL;
460 err = create_maps(Comm, B_file, B_row_map, B_col_map, B_range_map, B_domain_map);
461 if (err != 0) {
462 std::cout << "create_maps A returned err=="<<err<<std::endl;
463 return(err);
464 }
465
466 err = read_matrix(A_file, Comm, A_row_map, A_col_map,
467 A_range_map, A_domain_map, A);
468 if (err != 0) {
469 std::cout << "read_matrix A returned err=="<<err<<std::endl;
470 return(err);
471 }
472
473 err = read_matrix(B_file, Comm, B_row_map, B_col_map,
474 B_range_map, B_domain_map, B);
475 if (err != 0) {
476 std::cout << "read_matrix B returned err=="<<err<<std::endl;
477 return(-1);
478 }
479
480 const Epetra_Map* rowmap = transA ? &(A->DomainMap()) : &(A->RowMap());
481 const Epetra_Map* domainMap = transB ? &(B->RangeMap()) : &(B->DomainMap());
482 const Epetra_Map* rangeMap = transA ? &(A->DomainMap()) : &(A->RangeMap());
483
484
485 C = new Epetra_CrsMatrix(Copy, *rowmap, 1);
486
487 if(C->Comm().MyPID()) printf("transA = %d transB = %d\n",(int)transA,(int)transB);
488
489 err = EpetraExt::MatrixMatrix::Multiply(*A, transA, *B, transB, *C,do_FillComplete);
490 if (err != 0) {
491 std::cout << "err "<<err<<" from MatrixMatrix::Multiply"<<std::endl;
492 return(err);
493 }
494
495 if(!do_FillComplete) C->FillComplete(*domainMap,*rangeMap);
496
497
498// std::cout << "A: " << *A << std::endl << "B: "<<*B<<std::endl<<"C: "<<*C<<std::endl;
499 if (result_mtx_to_file) {
501 }
502
503 Epetra_Map* Cck_row_map = NULL;
504 Epetra_Map* Cck_col_map = NULL;
505 Epetra_Map* Cck_range_map = NULL;
506 Epetra_Map* Cck_domain_map = NULL;
507 err = create_maps(Comm, C_file, Cck_row_map, Cck_col_map,
508 Cck_range_map, Cck_domain_map);
509 if (err != 0) {
510 std::cout << "create_maps C returned err=="<<err<<std::endl;
511 return(err);
512 }
513
514 err = read_matrix(C_file, Comm, Cck_row_map, Cck_col_map,
515 Cck_range_map, Cck_domain_map, C_check);
516 if (err != 0) {
517 std::cout << "read_matrix C returned err=="<<err<<std::endl;
518 return(-1);
519 }
520
521 //now subtract our check-matrix from our result matrix (C) and
522 //the infinity-norm of that should be zero.
523 EpetraExt::MatrixMatrix::Add(*C, false, -1.0, *C_check, 1.0);
524
525 double inf_norm = C_check->NormInf();
526
527 int return_code = 0;
528
529 if (inf_norm < 1.e-13) {
530 if (localProc == 0 && verbose) {
531 std::cout << "Test Passed" << std::endl;
532 }
533 }
534 else {
535 return_code = -1;
536 if (localProc == 0) {
537 std::cout << "Test Failed ("<<filename<<"), inf_norm = " << inf_norm << std::endl;
538 }
539Comm.Barrier();
540std::cout << "C"<<std::endl;
541std::cout << *C<<std::endl;
542 }
543
544 delete A;
545 delete B;
546 delete C;
547 delete C_check;
548
549 delete A_row_map;
550 delete A_col_map;
551 delete A_range_map;
552 delete A_domain_map;
553
554 delete B_row_map;
555 delete B_col_map;
556 delete B_range_map;
557 delete B_domain_map;
558
559 delete Cck_row_map;
560 delete Cck_col_map;
561 delete Cck_range_map;
562 delete Cck_domain_map;
563
564 return(return_code);
565}
566
568 const std::string& input_file_name,
569 std::string& A_file,
570 bool& transA,
571 std::string& B_file,
572 bool& transB,
573 std::string& C_file)
574{
575 if (Comm.MyPID()==0) {
576 std::ifstream infile(input_file_name.c_str());
577 if (!infile) {
578 std::cout << "error opening input file " << input_file_name << std::endl;
579 return(-1);
580 }
581
582 const int linelen = 512;
583 char line[linelen];
584
585 infile.getline(line, linelen);
586 if (!infile.eof()) {
587 if (strchr(line, '#') == NULL) {
588 A_file = path+"/"+line;
589 }
590 }
591
592 infile.getline(line, linelen);
593 if (!infile.eof()) {
594 if (!strcmp(line, "TRANSPOSE")) {
595 transA = true;
596 }
597 else transA = false;
598 }
599
600 infile.getline(line, linelen);
601 if (!infile.eof()) {
602 if (strchr(line, '#') == NULL) {
603 B_file = path+"/"+line;
604 }
605 }
606
607 infile.getline(line, linelen);
608 if (!infile.eof()) {
609 if (!strcmp(line, "TRANSPOSE")) {
610 transB = true;
611 }
612 else transB = false;
613 }
614
615 infile.getline(line, linelen);
616 if (!infile.eof()) {
617 if (strchr(line, '#') == NULL) {
618 C_file = path+"/"+line;
619 }
620 }
621
622 broadcast_name(Comm, A_file);
623 broadcast_name(Comm, B_file);
624 broadcast_name(Comm, C_file);
625#ifdef EPETRA_MPI
626 int len = transA ? 1 : 0;
627 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
628 len = transB ? 1 : 0;
629 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
630#endif
631 }
632 else {
633 broadcast_name(Comm, A_file);
634 broadcast_name(Comm, B_file);
635 broadcast_name(Comm, C_file);
636#ifdef EPETRA_MPI
637 int len = 0;
638 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
639 transA = len==1 ? true : false;
640 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
641 transB = len==1 ? true : false;
642#endif
643 }
644
645 return(0);
646}
647
649 const std::string& input_file_name,
650 Epetra_Map*& row_map,
651 Epetra_Map*& col_map,
652 Epetra_Map*& range_map,
653 Epetra_Map*& domain_map)
654{
655 return( EpetraExt::MatrixMarketFileToBlockMaps64(input_file_name.c_str(),
656 Comm,
657 (Epetra_BlockMap*&)row_map,
658 (Epetra_BlockMap*&)col_map,
659 (Epetra_BlockMap*&)range_map,
660 (Epetra_BlockMap*&)domain_map)
661 );
662}
663
664int read_matrix(const std::string& filename,
665 Epetra_Comm& Comm,
666 const Epetra_Map* rowmap,
667 Epetra_Map* colmap,
668 const Epetra_Map* rangemap,
669 const Epetra_Map* domainmap,
670 Epetra_CrsMatrix*& mat)
671{
672 (void)Comm;
673 int err = EpetraExt::MatrixMarketFileToCrsMatrix64(filename.c_str(), *rowmap,
674 *rangemap, *domainmap, mat);
675
676 return(err);
677}
678
680 bool verbose)
681{
682 (void)verbose;
683 int thisproc = Comm.MyPID();
684 int numprocs = Comm.NumProc();
685
686 //only run this test on 2 procs
687 if (numprocs != 2) return(0);
688
689 //set up a row-std::map with 2 global elements,
690 //1 on each proc.
691 long long numGlobalRows = 2;
692 int numMyRows = 1;
693 long long myrow = 3;
694 if (thisproc == 1) myrow = 7;
695 Epetra_Map rowmap(numGlobalRows, numMyRows, &myrow, 0LL, Comm);
696
697 //set up a domain-map with columns 0 - 4 on proc 0,
698 //and columns 5 - 9 on proc 1.
699 long long numGlobalCols = 10;
700 int numMyCols = 5;
701 long long* mycols = new long long[numGlobalCols];
702 int i;
703 for(i=0; i<numGlobalCols; ++i) {
704 mycols[i] = i;
705 }
706
707 Epetra_Map domainmap(numGlobalCols, numMyCols, &(mycols[thisproc*numMyCols]),
708 0LL, Comm);
709
710 //now create matrices A, B and C with rowmap.
711 Epetra_CrsMatrix A(Copy, rowmap, 10);
712 Epetra_CrsMatrix B(Copy, rowmap, 10);
713 Epetra_CrsMatrix C(Copy, rowmap, 10);
714
715 double* coefs = new double[numGlobalCols];
716 for(i=0; i<numGlobalCols; ++i) {
717 coefs[i] = 1.0*i;
718 }
719
720 long long numCols = numGlobalCols - 2;
721 int offset = 0;
722 if (thisproc == 1) offset = 2;
723
724 int err = A.InsertGlobalValues(myrow, numCols, &coefs[offset], &mycols[offset]);
725
726 err += B.InsertGlobalValues(myrow, numMyCols, &(coefs[thisproc*numMyCols]),
727 &(mycols[thisproc*numMyCols]));
728
729 err += A.FillComplete(domainmap, rowmap);
730 err += B.FillComplete(domainmap, rowmap);
731
732 err += EpetraExt::MatrixMatrix::Multiply(A, false, B, true, C);
733
734 //std::cout << "two_proc_test, A: "<<std::endl;
735 //std::cout << A << std::endl;
736
737 //std::cout << "two_proc_test, B: "<<std::endl;
738 //std::cout << B << std::endl;
739
740 //std::cout << "two_proc_test, C: "<<std::endl;
741 //std::cout << C << std::endl;
742
743 if (C.NumGlobalNonzeros64() != 4) {
744 err += 1;
745 }
746
747 delete [] coefs;
748 delete [] mycols;
749
750 return(err);
751}
752
754{
755
756 const int magic_num = 3000;
757 // 2009/02/23: rabartl: If you are going to do a timing test you need to
758 // make this number adjustable form the command-line and you need to put in
759 // a real test that compares against hard numbers for pass/fail.
760
761 int localn = magic_num/Comm.NumProc();
762
764 Comm.MyPID(),
765 localn);
766
767 Epetra_CrsMatrix* C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
768
769 Epetra_Time timer(Comm);
770 double start_time = timer.WallTime();
771
772 int err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, false, *C);
773
774 int globaln = localn*Comm.NumProc();
775 if (verbose) {
776 std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A, time: "
777 << timer.WallTime()-start_time << std::endl;
778 }
779
780 C->FillComplete();
781
782 start_time = timer.WallTime();
783
784 err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, false, *C);
785
786 if (verbose) {
787 std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A, time: "
788 << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
789 }
790
791 delete C;
792
793 C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
794
795 start_time = timer.WallTime();
796
797 err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, true, *C);
798
799 if (verbose) {
800 std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A^T, time: "
801 << timer.WallTime()-start_time << std::endl;
802 }
803
804 C->FillComplete();
805
806 start_time = timer.WallTime();
807
808 err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, true, *C);
809
810 if (verbose) {
811 std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A^T, time: "
812 << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
813 }
814
815 delete C;
816
817 C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
818
819 start_time = timer.WallTime();
820
821 err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, false, *C);
822
823 if (verbose) {
824 std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A, time: "
825 << timer.WallTime()-start_time << std::endl;
826 }
827
828 C->FillComplete();
829
830 start_time = timer.WallTime();
831
832 err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, false, *C);
833
834 if (verbose) {
835 std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A, time: "
836 << timer.WallTime()-start_time << " (C already Filled)\n"<<std::endl;
837 }
838
839 delete C;
840
841 C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
842
843 start_time = timer.WallTime();
844
845 err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, true, *C);
846
847 if (verbose) {
848 std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A^T, time: "
849 << timer.WallTime()-start_time << std::endl;
850 }
851
852 C->FillComplete();
853
854 start_time = timer.WallTime();
855
856 err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, true, *C);
857
858 if (verbose) {
859 std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A^T, time: "
860 << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
861 }
862
863 delete C;
864
865 delete A;
866
867 return(err);
868}
869
871 int localProc,
872 int local_n,
873 bool callFillComplete,
874 bool symmetric)
875{
876 (void)localProc;
877#ifdef EPETRA_MPI
878 Epetra_MpiComm comm(MPI_COMM_WORLD);
879#else
881#endif
882 long long global_num_rows = numProcs*local_n;
883
884 Epetra_Map rowmap(global_num_rows, local_n, 0LL, comm);
885
886 int nnz_per_row = 9;
887 Epetra_CrsMatrix* matrix =
888 new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);
889
890 // Add rows one-at-a-time
891 double negOne = -1.0;
892 double posTwo = 2.0;
893 double val_L = symmetric ? negOne : -0.5;
894
895 for (int i=0; i<local_n; i++) {
896 long long GlobalRow = matrix->GRID64(i);
897 long long RowLess1 = GlobalRow - 1;
898 long long RowPlus1 = GlobalRow + 1;
899 long long RowLess5 = GlobalRow - 5;
900 long long RowPlus5 = GlobalRow + 5;
901 long long RowLess9 = GlobalRow - 9;
902 long long RowPlus9 = GlobalRow + 9;
903 long long RowLess24 = GlobalRow - 24;
904 long long RowPlus24 = GlobalRow + 24;
905 long long RowLess48 = GlobalRow - 48;
906 long long RowPlus48 = GlobalRow + 48;
907
908// if (!symmetric) RowLess5 -= 2;
909
910 if (RowLess48>=0) {
911 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess48);
912 }
913 if (RowLess24>=0) {
914 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess24);
915 }
916 if (RowLess9>=0) {
917 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess9);
918 }
919 if (RowLess5>=0) {
920 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess5);
921 }
922 if (RowLess1>=0) {
923 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess1);
924 }
925 if (RowPlus1<global_num_rows) {
926 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
927 }
928 if (RowPlus5<global_num_rows) {
929 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus5);
930 }
931 if (RowPlus9<global_num_rows) {
932 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus9);
933 }
934 if (RowPlus24<global_num_rows) {
935 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus24);
936 }
937 if (RowPlus48<global_num_rows) {
938 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus48);
939 }
940
941 matrix->InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
942 }
943
944 if (callFillComplete) {
945 int err = matrix->FillComplete();
946 if (err != 0) {
947 std::cout << "create_epetra_matrix: error in matrix.FillComplete()"
948 <<std::endl;
949 }
950 }
951
952 return(matrix);
953}
954
956{
957 int size = Comm.NumProc();
958 if (size != 2) return 0;
959
960 int rank = Comm.MyPID();
961
962 long long indexBase = 0;
963 long long numGlobalElements = 2;
964
965 Epetra_Map emap(numGlobalElements, indexBase, Comm);
966
967 Epetra_CrsMatrix A(Copy, emap, 0);
968
969 // 2x2 matrix:
970 // 3 4
971 // 1 2
972 std::vector<std::vector<double> > vals(numGlobalElements);
973 vals[0].push_back(3); vals[0].push_back(4);
974 vals[1].push_back(1); vals[1].push_back(2);
975
976 std::vector<long long> indices;
977 indices.push_back(0); indices.push_back(1);
978
979 for (int row=0; row<numGlobalElements; ++row) {
980 if ( A.MyGRID(row) )
981 A.InsertGlobalValues(row, numGlobalElements, &(vals[row][0]), &indices[0]);
982 }
983
984 A.FillComplete();
985
986 Epetra_CrsMatrix B(Copy, emap, 0);
987 EpetraExt::MatrixMatrix::Multiply(A, true, A, false, B);
988
989 // B = Transpose(A) x A should be
990 // 10 14
991 // 14 20
992 long long idx[2];
993 int tmp;
994 double val[2];
995
996 //for this little test, global_row == rank:
997 B.ExtractGlobalRowCopy(rank, 2, tmp, val, idx);
998
999 int test_result = 0;
1000
1001 if (rank == 0) {
1002 if (idx[0] == 0 && val[0] != 10.0) test_result = 1;
1003 if (idx[1] == 0 && val[1] != 10.0) test_result = 1;
1004 if (idx[0] == 1 && val[0] != 14.0) test_result = 1;
1005 if (idx[1] == 1 && val[1] != 14.0) test_result = 1;
1006 }
1007 else {
1008 if (idx[0] == 0 && val[0] != 14.0) test_result = 1;
1009 if (idx[1] == 0 && val[1] != 14.0) test_result = 1;
1010 if (idx[0] == 1 && val[0] != 20.0) test_result = 1;
1011 if (idx[1] == 1 && val[1] != 20.0) test_result = 1;
1012 }
1013
1014 int global_test_result = 0;
1015 Comm.SumAll(&test_result, &global_test_result, 1);
1016
1017 return global_test_result;
1018}
1019
Copy
static int Add(const Epetra_CrsMatrix &A, bool transposeA, double scalarA, Epetra_CrsMatrix &B, double scalarB)
Given Epetra_CrsMatrix objects A and B, form the sum B = a*A + b*B.
static int Multiply(const Epetra_CrsMatrix &A, bool transposeA, const Epetra_CrsMatrix &B, bool transposeB, Epetra_CrsMatrix &C, bool call_FillComplete_on_result=true, bool keep_all_hard_zeros=false)
Given Epetra_CrsMatrix objects A, B and C, form the product C = A*B.
int NumMyElements() const
virtual int NumProc() const=0
virtual int MyPID() const=0
virtual void Barrier() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
const Epetra_Map & RowMap() const
bool MyGRID(int GRID_in) const
int FillComplete(bool OptimizeDataStorage=true)
long long NumGlobalNonzeros64() const
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
const Epetra_Comm & Comm() const
const Epetra_Map & RangeMap() const
const Epetra_Map & DomainMap() const
long long GRID64(int LRID_in) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
double NormInf() const
int NumProc() const
int MyPID() const
double WallTime(void) const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int MatrixMarketFileToCrsMatrix64(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
int MatrixMarketFileToBlockMaps64(const char *filename, const Epetra_Comm &comm, Epetra_BlockMap *&rowmap, Epetra_BlockMap *&colmap, Epetra_BlockMap *&rangemap, Epetra_BlockMap *&domainmap)
Constructs row,col,range and domain maps from a matrix-market matrix file.
int RowMatrixToMatrixMarketFile(const char *filename, const Epetra_RowMatrix &A, const char *matrixName, const char *matrixDescription, bool writeHeader)
Writes an Epetra_RowMatrix object to a Matrix Market format file.
Epetra_Map * find_rows_containing_cols(const Epetra_CrsMatrix &M, const Epetra_Map &column_map)
int read_matrix_file_names(Epetra_Comm &Comm, const std::string &input_file_name, std::string &A_file, bool &transA, std::string &B_file, bool &transB, std::string &C_file)
int main(int argc, char **argv)
std::string path
int read_matrix(const std::string &filename, Epetra_Comm &Comm, const Epetra_Map *rowmap, Epetra_Map *colmap, const Epetra_Map *rangemap, const Epetra_Map *domainmap, Epetra_CrsMatrix *&mat)
int broadcast_name(Epetra_Comm &Comm, std::string &name)
int two_proc_test(Epetra_Comm &Comm, bool verbose=false)
int test_drumm1(Epetra_Comm &Comm)
int test_find_rows(Epetra_Comm &Comm)
int read_input_file(Epetra_Comm &Comm, const std::string &input_file_name, std::vector< std::string > &filenames)
Epetra_CrsMatrix * create_epetra_crsmatrix(int numProcs, int localProc, int local_n, bool callFillComplete=true, bool symmetric=true)
int expand_name_list(const char *newname, const char **&names, int &alloc_len, int &num_names)
int create_maps(Epetra_Comm &Comm, const std::string &input_file_name, Epetra_Map *&row_map, Epetra_Map *&col_map, Epetra_Map *&range_map, Epetra_Map *&domain_map)
int time_matrix_matrix_multiply(Epetra_Comm &Comm, bool verbose)
int run_test(Epetra_Comm &Comm, const std::string &filename, bool do_FillComplete, bool result_mtx_to_file=false, bool verbose=false)