[Trilinos-Users] matrix construction by Epetra

[Heroux, Mike]

If you want to continue down the path you are on, I suggest you use the std::map with keys that are the column entries and values that are the matrix coefficients.  The map can collect the values and then be unpacked to be insert.

This would be the declaration:
  std::map<int, double> curRowValues;

For each row, insert the values in the map in a loop.  This is the body:

      curRowValues[colGID] += value;

Make sure to clear the map in between rows.

Then to extract the values/indices, and insert them, you can do something like this:

  int i=0;
  std::map<int, double>::iterator pos;
  for (pos = curRowValues.begin(); pos != curRowValues.end(); ++pos) {
    indices[i] = pos->first;
    values[i++] = pos->second;
  }
   A->InsertGlobalValues(curRow, numEntries, values, indices);


Mike

From: <Heroux>, Michael A Heroux <mheroux at csbsju.edu<mailto:mheroux at csbsju.edu>>
Date: Thursday, January 16, 2014 7:33 AM
To: tin <tin at mma.cs.tsukuba.ac.jp<mailto:tin at mma.cs.tsukuba.ac.jp>>, "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: Re: [Trilinos-Users] matrix construction by Epetra

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>
┗━━━━━━━━━━━━━━━━━━━━━━━━━━━