[Trilinos-Users] matrix construction by Epetra

[Heroux, Mike]

I think the issue may be how the statement:

Z->SumIntoGlobalValues(MyGlobalElements[i], Num_B, Values_B, Indices_B);

Works.  If the entry you pass into this function is not already present in the matrix, the entry is ignored.  My guess is that you have entries in B that are not in A, so that these entries are not being kept.

I think your best option is to use the EpetraExt::MatrixMatrix::Add function.  It should work for you.

Also, to help debug, you should check the error codes from the functions you call.  They will tell you if errors are occurring.

I hope this helps.

Mike

From: tin <tin at mma.cs.tsukuba.ac.jp<mailto:tin at mma.cs.tsukuba.ac.jp>>
Reply-To: tin <tin at mma.cs.tsukuba.ac.jp<mailto:tin at mma.cs.tsukuba.ac.jp>>
Date: Thursday, January 16, 2014 12:22 AM
To: "trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>" <trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>>
Subject: [Trilinos-Users] matrix construction by Epetra


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,
Indices_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::endl;
std::cout << "Global Columns number: " << Z->NumGlobalCols() << std::endl;
std::cout << "Local Columns number: " << Z->NumMyCols() << std::endl;
}



┏┏┏┏━━━━━━━━━━━━━━━━━━━━━━━━
┏╋┏　　筑波大学大学院　　システム情報工学研究科　
┏┏■　　 CS専攻　 博士前期課程１年
┏　　　　　　情報数理研究室
┃　　　　　　　陳　玉蕊　（ちん　ぎょくずい）　　　　　
┃　　　　　　　e-mail: chen_yurui_cnu at yahoo.co.jp<mailto:chen_yurui_cnu at yahoo.co.jp>
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━