EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/MatrixMatrix/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 int numglobalrows = numprocs*numlocalrows;
252 Epetra_Map rowmap(numlocalrows*numprocs, 0, Comm);
253 Epetra_CrsMatrix matrix(Copy, rowmap, numglobalrows);
254
255 int err = 0;
256 int* cols = new int[numglobalrows];
257 double*vals = new double[numglobalrows];
258
259 for(int j=0; j<numglobalrows; ++j) {
260 cols[j] = j;
261 vals[j] = 1.0;
262 }
263
264 Epetra_Map colmap(-1, numglobalrows, cols, 0, Comm);
265
266 for(int i=0; i<numlocalrows; ++i) {
267 int row = 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 B 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 // Test Matrix Jacobi for non-transpose matrices
545 if(!transA && !transB && A->RowMap().SameAs(B->RowMap())) {
546 if(!Comm.MyPID()) std::cout<<"--Testing: Jacobi for same matrices"<<std::endl;
547 delete C; delete C_check;
548
549 double omega=1.0;
550 Epetra_Vector Dinv(B->RowMap());
551 Dinv.PutScalar(1.0);
552
553 // Jacobi version
554 C = new Epetra_CrsMatrix(Copy,B->RowMap(),0);
555 EpetraExt::MatrixMatrix::Jacobi(omega,Dinv,*A,*B,*C);
556
557 // Multiply + Add version
558 Dinv.PutScalar(omega);
559 Epetra_CrsMatrix * AB = new Epetra_CrsMatrix(Copy,B->RowMap(),0);
560 C_check = new Epetra_CrsMatrix(Copy,B->RowMap(),0);
561 EpetraExt::MatrixMatrix::Multiply(*A,false,*B,false,*AB);
562 AB->LeftScale(Dinv);
563 EpetraExt::MatrixMatrix::Add(*AB,false,-1.0,*B,false,1.0,C_check);
564
565 // Check the difference
566 EpetraExt::MatrixMatrix::Add(*C, false, -1.0, *C_check, 1.0);
567 C_check->FillComplete(B->DomainMap(),B->RangeMap());
568
569 // Error check
570 inf_norm = C_check->NormInf();
571 return_code = 0;
572 if (inf_norm < 1.e-13) {
573 if (localProc == 0 && verbose) {
574 std::cout << "Jacobi Test Passed" << std::endl;
575 }
576 }
577 else {
578 return_code = -1;
579 if (localProc == 0) {
580 std::cout << "Jacobi Test Failed ("<<filename<<"), inf_norm = " << inf_norm << std::endl;
581 }
582 Comm.Barrier();
583 std::cout << "C"<<std::endl;
584 std::cout << *C<<std::endl;
585 }
586
587 delete AB;
588 }
589
590
591 delete A;
592 delete B;
593 delete C;
594 delete C_check;
595
596 delete A_row_map;
597 delete A_col_map;
598 delete A_range_map;
599 delete A_domain_map;
600
601 delete B_row_map;
602 delete B_col_map;
603 delete B_range_map;
604 delete B_domain_map;
605
606 delete Cck_row_map;
607 delete Cck_col_map;
608 delete Cck_range_map;
609 delete Cck_domain_map;
610
611 return(return_code);
612}
613
615 const std::string& input_file_name,
616 std::string& A_file,
617 bool& transA,
618 std::string& B_file,
619 bool& transB,
620 std::string& C_file)
621{
622 if (Comm.MyPID()==0) {
623 std::ifstream infile(input_file_name.c_str());
624 if (!infile) {
625 std::cout << "error opening input file " << input_file_name << std::endl;
626 return(-1);
627 }
628
629 const int linelen = 512;
630 char line[linelen];
631
632 infile.getline(line, linelen);
633 if (!infile.eof()) {
634 if (strchr(line, '#') == NULL) {
635 A_file = path+"/"+line;
636 }
637 }
638
639 infile.getline(line, linelen);
640 if (!infile.eof()) {
641 if (!strcmp(line, "TRANSPOSE")) {
642 transA = true;
643 }
644 else transA = false;
645 }
646
647 infile.getline(line, linelen);
648 if (!infile.eof()) {
649 if (strchr(line, '#') == NULL) {
650 B_file = path+"/"+line;
651 }
652 }
653
654 infile.getline(line, linelen);
655 if (!infile.eof()) {
656 if (!strcmp(line, "TRANSPOSE")) {
657 transB = true;
658 }
659 else transB = false;
660 }
661
662 infile.getline(line, linelen);
663 if (!infile.eof()) {
664 if (strchr(line, '#') == NULL) {
665 C_file = path+"/"+line;
666 }
667 }
668
669 broadcast_name(Comm, A_file);
670 broadcast_name(Comm, B_file);
671 broadcast_name(Comm, C_file);
672#ifdef EPETRA_MPI
673 int len = transA ? 1 : 0;
674 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
675 len = transB ? 1 : 0;
676 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
677#endif
678 }
679 else {
680 broadcast_name(Comm, A_file);
681 broadcast_name(Comm, B_file);
682 broadcast_name(Comm, C_file);
683#ifdef EPETRA_MPI
684 int len = 0;
685 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
686 transA = len==1 ? true : false;
687 MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
688 transB = len==1 ? true : false;
689#endif
690 }
691
692 return(0);
693}
694
696 const std::string& input_file_name,
697 Epetra_Map*& row_map,
698 Epetra_Map*& col_map,
699 Epetra_Map*& range_map,
700 Epetra_Map*& domain_map)
701{
702 return( EpetraExt::MatrixMarketFileToBlockMaps(input_file_name.c_str(),
703 Comm,
704 (Epetra_BlockMap*&)row_map,
705 (Epetra_BlockMap*&)col_map,
706 (Epetra_BlockMap*&)range_map,
707 (Epetra_BlockMap*&)domain_map)
708 );
709}
710
711int read_matrix(const std::string& filename,
712 Epetra_Comm& Comm,
713 const Epetra_Map* rowmap,
714 Epetra_Map* colmap,
715 const Epetra_Map* rangemap,
716 const Epetra_Map* domainmap,
717 Epetra_CrsMatrix*& mat)
718{
719 (void)Comm;
720 int err = EpetraExt::MatrixMarketFileToCrsMatrix(filename.c_str(), *rowmap,
721 *rangemap, *domainmap, mat);
722
723 return(err);
724}
725
727 bool verbose)
728{
729 (void)verbose;
730 int thisproc = Comm.MyPID();
731 int numprocs = Comm.NumProc();
732
733 //only run this test on 2 procs
734 if (numprocs != 2) return(0);
735
736 //set up a row-std::map with 2 global elements,
737 //1 on each proc.
738 int numGlobalRows = 2;
739 int numMyRows = 1;
740 int myrow = 3;
741 if (thisproc == 1) myrow = 7;
742 Epetra_Map rowmap(numGlobalRows, numMyRows, &myrow, 0, Comm);
743
744 //set up a domain-map with columns 0 - 4 on proc 0,
745 //and columns 5 - 9 on proc 1.
746 int numGlobalCols = 10;
747 int numMyCols = 5;
748 int* mycols = new int[numGlobalCols];
749 int i;
750 for(i=0; i<numGlobalCols; ++i) {
751 mycols[i] = i;
752 }
753
754 Epetra_Map domainmap(numGlobalCols, numMyCols, &(mycols[thisproc*numMyCols]),
755 0, Comm);
756
757 //now create matrices A, B and C with rowmap.
758 Epetra_CrsMatrix A(Copy, rowmap, 10);
759 Epetra_CrsMatrix B(Copy, rowmap, 10);
760 Epetra_CrsMatrix C(Copy, rowmap, 10);
761
762 double* coefs = new double[numGlobalCols];
763 for(i=0; i<numGlobalCols; ++i) {
764 coefs[i] = 1.0*i;
765 }
766
767 int numCols = numGlobalCols - 2;
768 int offset = 0;
769 if (thisproc == 1) offset = 2;
770
771 int err = A.InsertGlobalValues(myrow, numCols, &coefs[offset], &mycols[offset]);
772
773 err += B.InsertGlobalValues(myrow, numMyCols, &(coefs[thisproc*numMyCols]),
774 &(mycols[thisproc*numMyCols]));
775
776 err += A.FillComplete(domainmap, rowmap);
777 err += B.FillComplete(domainmap, rowmap);
778
779 err += EpetraExt::MatrixMatrix::Multiply(A, false, B, true, C);
780
781 //std::cout << "two_proc_test, A: "<<std::endl;
782 //std::cout << A << std::endl;
783
784 //std::cout << "two_proc_test, B: "<<std::endl;
785 //std::cout << B << std::endl;
786
787 //std::cout << "two_proc_test, C: "<<std::endl;
788 //std::cout << C << std::endl;
789
790 if (C.NumGlobalNonzeros() != 4) {
791 err += 1;
792 }
793
794 delete [] coefs;
795 delete [] mycols;
796
797 return(err);
798}
799
801{
802
803 const int magic_num = 3000;
804 // 2009/02/23: rabartl: If you are going to do a timing test you need to
805 // make this number adjustable form the command-line and you need to put in
806 // a real test that compares against hard numbers for pass/fail.
807
808 int localn = magic_num/Comm.NumProc();
809
811 Comm.MyPID(),
812 localn);
813
814 Epetra_CrsMatrix* C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
815
816 Epetra_Time timer(Comm);
817 double start_time = timer.WallTime();
818
819 int err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, false, *C);
820
821 int globaln = localn*Comm.NumProc();
822 if (verbose) {
823 std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A, time: "
824 << timer.WallTime()-start_time << std::endl;
825 }
826
827 C->FillComplete();
828
829 start_time = timer.WallTime();
830
831 err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, false, *C);
832
833 if (verbose) {
834 std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A, time: "
835 << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
836 }
837
838 delete C;
839
840 C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
841
842 start_time = timer.WallTime();
843
844 err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, true, *C);
845
846 if (verbose) {
847 std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A^T, time: "
848 << timer.WallTime()-start_time << std::endl;
849 }
850
851 C->FillComplete();
852
853 start_time = timer.WallTime();
854
855 err = EpetraExt::MatrixMatrix::Multiply(*A, false, *A, true, *C);
856
857 if (verbose) {
858 std::cout << "size: " << globaln << "x"<<globaln<<", C = A*A^T, time: "
859 << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
860 }
861
862 delete C;
863
864 C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
865
866 start_time = timer.WallTime();
867
868 err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, false, *C);
869
870 if (verbose) {
871 std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A, time: "
872 << timer.WallTime()-start_time << std::endl;
873 }
874
875 C->FillComplete();
876
877 start_time = timer.WallTime();
878
879 err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, false, *C);
880
881 if (verbose) {
882 std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A, time: "
883 << timer.WallTime()-start_time << " (C already Filled)\n"<<std::endl;
884 }
885
886 delete C;
887
888 C = new Epetra_CrsMatrix(Copy, A->RowMap(), 0);
889
890 start_time = timer.WallTime();
891
892 err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, true, *C);
893
894 if (verbose) {
895 std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A^T, time: "
896 << timer.WallTime()-start_time << std::endl;
897 }
898
899 C->FillComplete();
900
901 start_time = timer.WallTime();
902
903 err = EpetraExt::MatrixMatrix::Multiply(*A, true, *A, true, *C);
904
905 if (verbose) {
906 std::cout << "size: " << globaln << "x"<<globaln<<", C = A^T*A^T, time: "
907 << timer.WallTime()-start_time << " (C already Filled)\n" <<std::endl;
908 }
909
910 delete C;
911
912 delete A;
913
914 return(err);
915}
916
918 int localProc,
919 int local_n,
920 bool callFillComplete,
921 bool symmetric)
922{
923 (void)localProc;
924#ifdef EPETRA_MPI
925 Epetra_MpiComm comm(MPI_COMM_WORLD);
926#else
928#endif
929 int global_num_rows = numProcs*local_n;
930
931 Epetra_Map rowmap(global_num_rows, local_n, 0, comm);
932
933 int nnz_per_row = 9;
934 Epetra_CrsMatrix* matrix =
935 new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);
936
937 // Add rows one-at-a-time
938 double negOne = -1.0;
939 double posTwo = 2.0;
940 double val_L = symmetric ? negOne : -0.5;
941
942 for (int i=0; i<local_n; i++) {
943 int GlobalRow = matrix->GRID(i);
944 int RowLess1 = GlobalRow - 1;
945 int RowPlus1 = GlobalRow + 1;
946 int RowLess5 = GlobalRow - 5;
947 int RowPlus5 = GlobalRow + 5;
948 int RowLess9 = GlobalRow - 9;
949 int RowPlus9 = GlobalRow + 9;
950 int RowLess24 = GlobalRow - 24;
951 int RowPlus24 = GlobalRow + 24;
952 int RowLess48 = GlobalRow - 48;
953 int RowPlus48 = GlobalRow + 48;
954
955// if (!symmetric) RowLess5 -= 2;
956
957 if (RowLess48>=0) {
958 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess48);
959 }
960 if (RowLess24>=0) {
961 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess24);
962 }
963 if (RowLess9>=0) {
964 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess9);
965 }
966 if (RowLess5>=0) {
967 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess5);
968 }
969 if (RowLess1>=0) {
970 matrix->InsertGlobalValues(GlobalRow, 1, &val_L, &RowLess1);
971 }
972 if (RowPlus1<global_num_rows) {
973 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
974 }
975 if (RowPlus5<global_num_rows) {
976 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus5);
977 }
978 if (RowPlus9<global_num_rows) {
979 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus9);
980 }
981 if (RowPlus24<global_num_rows) {
982 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus24);
983 }
984 if (RowPlus48<global_num_rows) {
985 matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus48);
986 }
987
988 matrix->InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
989 }
990
991 if (callFillComplete) {
992 int err = matrix->FillComplete();
993 if (err != 0) {
994 std::cout << "create_epetra_matrix: error in matrix.FillComplete()"
995 <<std::endl;
996 }
997 }
998
999 return(matrix);
1000}
1001
1003{
1004 int size = Comm.NumProc();
1005 if (size != 2) return 0;
1006
1007 int rank = Comm.MyPID();
1008
1009 int indexBase = 0;
1010 int numGlobalElements = 2;
1011
1012 Epetra_Map emap(numGlobalElements, indexBase, Comm);
1013
1014 Epetra_CrsMatrix A(Copy, emap, 0);
1015
1016 // 2x2 matrix:
1017 // 3 4
1018 // 1 2
1019 std::vector<std::vector<double> > vals(numGlobalElements);
1020 vals[0].push_back(3); vals[0].push_back(4);
1021 vals[1].push_back(1); vals[1].push_back(2);
1022
1023 std::vector<int> indices;
1024 indices.push_back(0); indices.push_back(1);
1025
1026 for (int row=0; row<numGlobalElements; ++row) {
1027 if ( A.MyGRID(row) )
1028 A.InsertGlobalValues(row, numGlobalElements, &(vals[row][0]), &indices[0]);
1029 }
1030
1031 A.FillComplete();
1032
1033 Epetra_CrsMatrix B(Copy, emap, 0);
1034 EpetraExt::MatrixMatrix::Multiply(A, true, A, false, B);
1035
1036 // B = Transpose(A) x A should be
1037 // 10 14
1038 // 14 20
1039 int idx[2];
1040 int tmp;
1041 double val[2];
1042
1043 //for this little test, global_row == rank:
1044 B.ExtractGlobalRowCopy(rank, 2, tmp, val, idx);
1045
1046 int test_result = 0;
1047
1048 if (rank == 0) {
1049 if (idx[0] == 0 && val[0] != 10.0) test_result = 1;
1050 if (idx[1] == 0 && val[1] != 10.0) test_result = 1;
1051 if (idx[0] == 1 && val[0] != 14.0) test_result = 1;
1052 if (idx[1] == 1 && val[1] != 14.0) test_result = 1;
1053 }
1054 else {
1055 if (idx[0] == 0 && val[0] != 14.0) test_result = 1;
1056 if (idx[1] == 0 && val[1] != 14.0) test_result = 1;
1057 if (idx[0] == 1 && val[0] != 20.0) test_result = 1;
1058 if (idx[1] == 1 && val[1] != 20.0) test_result = 1;
1059 }
1060
1061 int global_test_result = 0;
1062 Comm.SumAll(&test_result, &global_test_result, 1);
1063
1064 return global_test_result;
1065}
1066
Copy
static int Jacobi(double omega, const Epetra_Vector &Dinv, const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B, Epetra_CrsMatrix &C, bool call_FillComplete_on_result=true)
Given Epetra_CrsMatrix objects A, B and C, and Epetra_Vector Dinv, form the product C = (I-omega * Di...
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.
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)
int NumGlobalNonzeros() 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
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int NumProc() const
int MyPID() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int MatrixMarketFileToBlockMaps(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 MatrixMarketFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
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)