[Trilinos-Users] problems with transpose, inverse and products with Epetra_CRSMatrix

Fri Sep 12 08:54:18 MDT 2014

Small error below

  sum +=curRowVals[i];

Should be
  sum += curRowVals[j]*curRowVals[j];

Please note that I did not compile or test this code, so be suspicious of
details :).

Mike

On 9/12/14, 9:49 AM, "Heroux, Mike" <MHeroux at CSBSJU.EDU> wrote:

>Antonella,
>
>First, although I don’t think you will need it for this situation, you can
>enable EpetraExt using Cmake ENABLE option and then re-doing the make
>install.
>
>Second, forming B can be done much more efficiently than computing the
>entire normal matrix.  You really just need the diagonal of the matrix,
>which is obtained entry by entry from the sum of squares of the matrix
>rows.  The following code should be a good start (although I have not
>tested it) for computing diag(A^T A):
>
>Assuming A is your Epetra_CrsMatrix object and diagATA is an Epetra_Vector
>defined using the RowMap of A (Epetra_Vector diagATA(A.RowMap()))
>
> this function should compute diag(A^T A):
>
>  void computeB(const Epetra_CrsMatrix & A, Epetra_Vector & diagATA) {
>
>// Preconditions: RowMap of A == DomainMap of A == RangeMap of A (These
>conditions are typically true)
>
>    for(int i=0; i<A.NumMyRows(); i++) {
>	int numEntries = A.Graph().NumMyIndices(i);
>	double * curRowVals = A.Values(i);
>	double sum = 0.0;
>	for(int j=0; j< numEntries; j++)
>		sum += curRowVals[j];
>	diagATA[i] = sum;
>    }
>  return;
>  }
>
>On exit, diagATA will be a vector containing the diagonal elements of A^T
>A.
>
>
>Next you can compute B^T = A * diag(A^T A) as follows:
>
>  Epetra_CrsMatrix BT(A); // makes copy of A into B^T
>  BT.RightScale(diagATA); // Scale B in-place with diag(A^T A)
>
>At this point you have the elements you need for working with the system
>BAx=Bb.  
>
>You can continue by explicitly forming BA and Bb, if you like (not
>necessarily recommended, but might be OK).
>
>BA is formed by calling the MatrixMatrix functions Alicia recommended from
>EpetraExt.
>
>Bb is formed by:
>
>  Epetra_Vector Bb(BT.DomainMap());
>  BT.Multiply(true, b, Bb);
>
>Hopefully this is enough information to get you started.
>
>Mike
>
>On 9/11/14, 7:10 AM, "Antonella Longo" <longo at pi.ingv.it> wrote:
>
>>Thank you for your prompt answer.
>>In this email we try to explain the complete problem that we must solve.
>>
>>We have a non-square system of equation for which we found a
>>preconditioner
>>to make it square so that it can be solved with GMRES.
>>
>>Our initial non-square system is Ax=b.
>>The preconditioner requires the computation of the matrix B
>>that is $B=(diag(A^T A))^{-1}A^T$.
>>A is an Epetra_CRSMatrix and we have its map.
>>The final system to be solved will be BAx=Bb.
>>
>>We were not successful to make this B matrix defined above.
>>
>>Another problem that we have is that even though the EpetraExt package
>>is installed, and the proper ".h" files are included, the compiler does
>>not 
>>recognize the functions that you suggested us.
>>
>>We are wondering if you could suggest us the right functions
>>to solve the above problem.
>>
>>-----------------------------------------------
>>Antonella Longo
>>Istituto Nazionale di Geofisica e Vulcanologia
>>Sezione di Pisa
>>Via Uguccione della Faggiola, 32
>>Pisa, Italy
>>Tel. office: 050 8311939
>>Cell.: 339 7532089
>>
>>----- Original Message -----
>>From: "Mike Heroux" <MHeroux at CSBSJU.EDU>
>>To: "Antonella Longo" <longo at pi.ingv.it>, "Deepak Garg"
>><deepak.garg at pi.ingv.it>
>>Cc: trilinos-users at software.sandia.gov, "Alicia Marie Klinvex"
>><aklinvex at purdue.edu>
>>Sent: Wednesday, September 10, 2014 11:17:47 PM
>>Subject: Re: [Trilinos-Users] problems with transpose, inverse and
>>products with Epetra_CRSMatrix
>>
>>In addition to the pointers that Alicia provided, you can compute the
>>matrix-transpose times vector product using the specific Epetra_CrsMatrix
>>method called Multiply().  This method has a bool input argument that
>>indicates whether or not to apply the matrix-transpose operation:
>>
>>http://trilinos.org/docs/dev/packages/epetra/doc/html/classEpetra__CrsMat
>>r
>>i
>>x.html#ae97a6c65b11aa29b3bb25d9788c56b3a
>>
>>In general, it is unusual (but not necessarily wrong) to explicitly form
>>the transpose of a matrix.  In many situations, just being able to apply
>>the action of the matrix transpose is sufficient.
>>
>>Following up on Alicia¹s comments on inverses, it is almost never best to
>>explicitly compute the inverse of a sparse matrix, unless it is small or
>>has very special structure.
>>
>>Perhaps if you describe the algorithms you are trying to implement we can
>>suggest some other options.
>>
>>Mike
>>
>>On 9/10/14, 8:10 AM, "Alicia Marie Klinvex" <aklinvex at purdue.edu> wrote:
>>
>>>Hello Antonella and Deepak,
>>>
>>>The matrix multiplication is in the EpetraExt package:
>>>http://trilinos.org/docs/dev/packages/epetraext/doc/html/classEpetraExt_
>>>1
>>>_
>>>1MatrixMatrix.html
>>>
>>>I believe the transpose function you seek is also in EpetraExt:
>>>http://trilinos.org/docs/dev/packages/epetraext/doc/html/classEpetraExt_
>>>1
>>>_
>>>1RowMatrix__Transpose.html
>>>(I could be wrong about that one; I've never used that particular
>>>function before.)
>>>
>>>The matrix vector multiplication is in Epetra (called Apply):
>>>http://trilinos.org/docs/dev/packages/epetra/doc/html/classEpetra__CrsMa
>>>t
>>>r
>>>ix.html#a77c4232f1b3d32ff76504ed314df5f76
>>>
>>>Do you absolutely need to compute the inverse of your matrix, or do you
>>>just want to apply the inverse to a vector?  I'm afraid I don't know
>>>whether Trilinos even has a function to invert a sparse matrix, but if
>>>you were only interested in solving linear systems, it certainly has
>>>functions for that.
>>>
>>>Best wishes,
>>>Alicia Klinvex
>>>
>>>
>>>
>>>----- Original Message -----
>>>From: "Antonella Longo" <longo at pi.ingv.it>
>>>To: trilinos-users at software.sandia.gov
>>>Cc: "Deepak Garg" <deepak.garg at pi.ingv.it>
>>>Sent: Wednesday, September 10, 2014 9:10:46 AM
>>>Subject: [Trilinos-Users] problems with transpose,	inverse and products
>>>with Epetra_CRSMatrix
>>>
>>>We are making computations with Epetra_CRSMatrices,
>>>but we could not find any mathematical functions for
>>>these matrices.
>>>In particular we need to compute the transpose and
>>>the inverse of an Epetra_CRSMatrix, to calculate
>>>the product of two Epetra_CRSMatrices, and the
>>>product of an Epetra_CRSMatrix with an Epetra_Vector.
>>>Does there exist any standard functions for these operations?
>>>We checked the functions of the above mentioned classes but we found
>>>nothing. 
>>>
>>>Thanks for your kind attention.
>>>
>>>Antonella Longo and Deepak Garg
>>>
>>>
>>>
>>>-----------------------------------------------
>>>Antonella Longo
>>>Istituto Nazionale di Geofisica e Vulcanologia
>>>Sezione di Pisa
>>>Via Uguccione della Faggiola, 32
>>>Pisa, Italy
>>>Tel. office: 050 8311939
>>>Cell.: 339 7532089
>>>
>>>_______________________________________________
>>>Trilinos-Users mailing list
>>>Trilinos-Users at software.sandia.gov
>>>https://software.sandia.gov/mailman/listinfo/trilinos-users
>>>_______________________________________________
>>>Trilinos-Users mailing list
>>>Trilinos-Users at software.sandia.gov
>>>https://software.sandia.gov/mailman/listinfo/trilinos-users
>>
>
>
>_______________________________________________
>Trilinos-Users mailing list
>Trilinos-Users at software.sandia.gov
>https://software.sandia.gov/mailman/listinfo/trilinos-users