[Trilinos-Users] Problem with EpetraExt::MatrixMatrix::Multiply

[Williams, Alan B]

Maurizio,
Sorry about the delay in looking into this issue.
I reproduced your problem, here are a couple of points:


1.       The MatrixMatrix::Multiply function forms C = A*A,  *not* C = A*A + C.

2.       The MatrixMatrix::Multiply function returns an error-code which is 0 if successful, nonzero if some error occurred. In your case, it was returning a nonzero error-code. The issue is that your C matrix has FillComplete() called before being passed to the Multiply function, and so the Multiply function was unable to add additional nonzero elements to C. Note that C needs to contain more nonzeros than A in this case.

3.       If you call FillComplete() on A, create C using the copy-constructor to copy A, then C has FillComplete() called also.

4.       In your code, create Matrix_C like this:
Epetra_CrsMatrix Matrix_C(Copy,Map_func, 0);
then call the MatrixMatrix::Multiply function, and Matrix_C will contain the correct results.

5.       Check the return-code from Multiply to be sure that it returns 0 indicating success. Admittedly, the Multiply function is more fragile than it should be. A parallel matrix-matrix multiply that handles rectangular matrices is somewhat complicated.



Alan



From: trilinos-users-bounces at software.sandia.gov [mailto:trilinos-users-bounces at software.sandia.gov] On Behalf Of bender calculon
Sent: Sunday, September 20, 2009 4:32 AM
To: trilinos-users at software.sandia.gov
Subject: [Trilinos-Users] Problem with EpetraExt::MatrixMatrix::Multiply

Hi,

I have this problem with EpetraExt:MatrixMatrix::Multiply. I want to multiply two Epetra_CrsMatrix get from Trilinos_Util::CrsMatrixGallery.
If I multiply two full matrix ( for example Hilber Matrix, Cauchy matrix or kms matrix ) or two diagonal matrix I obtain a correct result, but if I multiply two sparse matrix i obtain an absurd result.

I explain my problem with a part of my code:

Trilinos_Util::CrsMatrixGallery Problem(problem, Comm);   // problem = "laplace_2d"
            Problem.Set("problem_size",size);                               // size = 4
            Problem.Set("map_type","linear");
            read_A = Problem.GetMatrix();                                   // Epetra_CrsMatrix *read_A
            Num_Glob = read_A->Map().NumGlobalElements();    // int Num_Glob
            Epetra_Map Map_func(Num_Glob,0,Comm);
            Epetra_Export exporter(read_A -> Map(), Map_func);     // I want to export  data from *read_A to matrix_A
            Epetra_CrsMatrix Matrix_A(Copy, Map_func, 0);
                        Matrix_A.Export(*read_A,exporter, Zero);
                        Matrix_A.FillComplete();
            cout << Matrix_A;
            Epetra_CrsMatrix Matrix_C(Matrix_A);                             // copy of Matrix_A;
            Matrix_C.PutScalar(0);
            Matrix_C.FillComplete();
            EpetraExt::MatrixMatrix::Multiply(Matrix_A,false,Matrix_A,false,Matrix_C);   // This operation make Matrix_C = Matrix_A*Matrix_A +C not C = A*A

If i put problem = "laplace_2d" and "problem_size" = 4 I obtain this result:

./ex2.exe -problem laplace_2d -size 4
Epetra::CrsMatrix

Number of Global Rows        = 4
Number of Global Cols        = 4
Number of Global Diagonals   = 4
Number of Global Nonzeros    = 12
Global Maximum Num Entries   = 3

Number of My Rows        = 4
Number of My Cols        = 4
Number of My Diagonals   = 4
Number of My Nonzeros    = 12
My Maximum Num Entries   = 3

   Processor    Row Index    Col Index           Value
       0             0             0                      17
       0             0             1                      -8
       0             0             2                      -4
       0             1             0                       0
       0             1             1                       0
       0             1             3                       0
       0             2             0                       0
       0             2             2                       0
       0             2             3                       0
       0             3             1                       0
       0             3             2                       0
       0             3             3                       0

and you can see that is not correct.

I have trilinos release 9.0.

T

Best regards

Maurizio
________________________________
Windows Live Messenger Windows Live Messenger, l'unico e l'originale. Diffida dalle imitazioni!<http://www.messenger.it>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20090924/4e0550f6/attachment.html