[Trilinos-Users] matrix construction by Epetra

tin tin at mma.cs.tsukuba.ac.jp
Thu Jan 16 01:22:58 MST 2014

Hi everyone, I'm new to Trilinos and Epetra and now working on implementation of one
algorithm. Before computation, I have to read matrices from .rsa files
and then do some easy matrix computaion to get a new matrix for next
computation. The problem appears in the construction of this new matrix. I want to compute the new matrix by : Z = wB - A, where A and B are
imput matrices and w is parameter. I've tried a lot of ways, but cannot
get correct answer. I think that the problem may in Map construction,
but cannot make sure. Here is part of my program. Thank you very much!
By the way, I also tried EpetraExt::MatrixMatrix:Add, but still stopped in FillComplete(). typedef double ST;
typedef Teuchos::ScalarTraits<ST> SCT;
typedef SCT::magnitudeType MT;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
typedef Belos::MultiVecTraits<ST,MV> MVT;
typedef Belos::OperatorTraits<ST,MV,OP> OPT; using Teuchos::ParameterList;
using Teuchos::RCP;
using Teuchos::rcp; bool proc_verbose = false;
if (MyPID == 0) proc_verbose = true; Trilinos_Util::CrsMatrixGallery A0("hb", Comm);
A0.Set("matrix_name", "bcsstk08.rsa");
A0.Set("map_type", "linear");
Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp(A0.GetMatrix(), false);
Trilinos_Util::CrsMatrixGallery B0("hb", Comm);
B0.Set("matrix_name", "bcsstm08.rsa");
B0.Set("map_type", "greedy");
Teuchos::RCP<Epetra_CrsMatrix> B = Teuchos::rcp(B0.GetMatrix(), false); // Z = wB - A
RCP<Epetra_Map> Map = rcp( new Epetra_Map(A->RowMap()) );
int NumMyElements = Map->NumMyElements();
int *MyGlobalElements = Map->MyGlobalElements();
const int Length = A->NumGlobalRows();
int *NumNz = new int[NumMyElements];
for (int i = 0; i < NumMyElements; i++)
NumNz[i] = A->NumMyEntries(i);
RCP<Epetra_CrsMatrix> Z = rcp( new Epetra_CrsMatrix(Copy, *Map, NumNz) ); double theta;
MT omega, pp;
double center = 1.0, radius = 0.1;
int N = 32;
theta = ((2*M_PI)/N)*(1+0.5);
pp = cos(theta);
omega = center + radius*pp;
double *Values_A = new double[Length];
double *Values_B = new double[Length];
int *Indices_A = new int[Length];
int *Indices_B = new int[Length];
int Num_A = 0, Num_B = 0; for (int i = 0; i < NumMyElements; i++)
A->ExtractGlobalRowCopy(MyGlobalElements[i], Length, Num_A, Values_A,
B->ExtractGlobalRowCopy(MyGlobalElements[i], Length, Num_B, Values_B,
Indices_B); for (int j = 0; j < Num_A; j++)
Values_A[j] = -1.0 * Values_A[j]; Z->InsertGlobalValues(MyGlobalElements[i], Num_A, Values_A, Indices_A);
if (Num_B == 0) continue;
for (int j = 0; j < Num_B; j++)
Values_B[j] = omega * Values_B[j]; Z->SumIntoGlobalValues(MyGlobalElements[i], Num_B, Values_B, Indices_B);
delete Values_A; delete Indices_A;
delete Values_B; delete Indices_B;
delete NumNz;
Z->FillComplete(); if (proc_verbose)
std::cout << "One-norm: " << Z->NormOne() << std::endl;
std::cout << "Global Nonzeros number: " << Z->NumGlobalNonzeros() <<
std::cout << "Global Columns number: " << Z->NumGlobalCols() << std::endl;
std::cout << "Local Columns number: " << Z->NumMyCols() << std::endl;
┏╋┏  筑波大学大学院  システム情報工学研究科 
┏┏■       CS専攻  博士前期課程1年
┏      情報数理研究室
┃       陳 玉蕊 (ちん ぎょくずい)     
┃       e-mail: chen_yurui_cnu at yahoo.co.jp
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20140116/7cb39e09/attachment.html 

More information about the Trilinos-Users mailing list