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

Fri Sep 12 09:52:04 MDT 2014

Sorry. Yes, of course.  But this is still possible.   You just need to
compute the reciprocal values on BT:

 BT.Reciprocal(BT);

Mike

On 9/12/14, 10:17 AM, "Deepak Garg" <deepak.garg at pi.ingv.it> wrote:

>Dear Mike
>
>Thank you very much for your kind reply.
>
>I have a little doubt. We need inverse of diag(A^T A).
>
>B^T should be 
>
>B^T = A * (diag(A^T A))^{-1}
>
>How will we compute this inverse?
>
>Thank you 
>Deepak
>
>----- Original Message -----
>From: "Mike Heroux" <MHeroux at CSBSJU.EDU>
>To: "Mike Heroux" <MHeroux at CSBSJU.EDU>, "Antonella Longo"
><longo at pi.ingv.it>
>Cc: trilinos-users at software.sandia.gov, "Deepak Garg"
><deepak.garg at pi.ingv.it>
>Sent: Friday, September 12, 2014 4:54:18 PM
>Subject: Re: [Trilinos-Users] problems with transpose, inverse and
>products with Epetra_CRSMatrix
>
>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__CrsMa
>>>t
>>>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__CrsM
>>>>a
>>>>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