[Trilinos-Users] Re: multiply two Epetra_CrsMatrix and get back the diagonal

[Kurt Stokbro]

Hi Michael,
One of the matrixes is symmetric, so I have implemented the attached
routine, that works in seriel. I was wondering if there is a routine
that works in parallel and that is more optimal.

Regards,
Kurt





On 10/30/07, Michael A Heroux <maherou at sandia.gov> wrote:
> Kurt,
>
> Given the row bias of Epetra_CrsMatrix, going down the columns of b is a bit
> tricky.  There are two basic approaches:
>
> - Transpose b all at once.  This would be most time-efficient, but take
> extra memory at least temporarily.
> - Grab a column of b at a time.  This would be most space efficient, but
> take more time.
>
> Do you have a sense of which is most acceptable?
>
> Mike
>
>
> On 10/29/07 2:19 PM, "Kurt Stokbro" <kurt.stokbro at gmail.com> wrote:
>
> > I need a subroutine that can
> > Epetra_Vector multiplyDiag(Epetra_CrsMatrix& a, Epetra_CrsMatrix& b)
> >  {
> >    for(int j=0;j<a.numRows;j++)
> >    {
> >          v(j)=0.0
> >         for(int i=0;i<a.numCols;i++)
> >                  v(j)+=a(j,i)*b(i,j)
> >   }
> >    return v
> > }
> >
> > Does anybody have an efficient routine for this?
> > Thanks,
> > Kurt Stokbro
> >
> > _______________________________________________
> > Trilinos-Users mailing list
> > Trilinos-Users at software.sandia.gov
> > http://software.sandia.gov/mailman/listinfo/trilinos-users
>
>
>
-------------- next part --------------
typedef Teuchos::SerialDenseVector<int,double> TVector;


  TVector multiplyDiag( Epetra_CrsMatrix & G, Epetra_CrsMatrix & S)
  {
    int numOrb= G.NumGlobalRows();
    int MaxNumIndicesG = G.MaxNumEntries();
    int extractedNumEntriesG;
    int * extractedIndicesG  = new int[MaxNumIndicesG];
    double * extractedValuesG  = new double[MaxNumIndicesG];
    int MaxNumIndicesS = S.MaxNumEntries();
    int extractedNumEntriesS;
    int * extractedIndicesS  = new int[MaxNumIndicesS];
    double * extractedValuesS  = new double[MaxNumIndicesS];

    TVector diag(numOrb) ;
    if(G.NumGlobalNonzeros()==S.NumGlobalNonzeros())
      {
	for(int i=0;i<numOrb;i++)
	  {
	    int err = G.ExtractGlobalRowCopy(
		    i,                     // input
                    MaxNumIndicesG,         //
                    extractedNumEntriesG,   // output
                    extractedValuesG,       // output
                    extractedIndicesG);     // output
	assert(err == 0);
	err = S.ExtractGlobalRowCopy(
		    i,                     // input
                    MaxNumIndicesS,         //
                    extractedNumEntriesS,   // output
                    extractedValuesS,       // output
                    extractedIndicesS);     // output
	assert(err == 0);


	for(int j= 0; j < extractedNumEntriesG; j++)
	  diag(i)+=extractedValuesG[j]*extractedValuesS[j];
	  }
      }
    else
      {
	for(int i=0;i<numOrb;i++)
	  {
	    
	    int err = G.ExtractGlobalRowCopy(
		    i,                     // input
                    MaxNumIndicesG,         //
                    extractedNumEntriesG,   // output
                    extractedValuesG,       // output
                    extractedIndicesG);     // output
	assert(err == 0);
	err = S.ExtractGlobalRowCopy(
		    i,                     // input
                    MaxNumIndicesS,         //
                    extractedNumEntriesS,   // output
                    extractedValuesS,       // output
                    extractedIndicesS);     // output
	assert(err == 0);
	int jS=0;
	for(int jG = 0; jG < extractedNumEntriesG; jG++)
            {
                int jiG = extractedIndicesG[jG];
		while(extractedIndicesS[jS]<jiG && jS<extractedNumEntriesS-1)
		  jS++;
		if (extractedIndicesS[jS]==jiG )
		  diag(i)+=extractedValuesG[jG]*extractedValuesS[jS];
	    }
	  }
      }

    delete [] extractedIndicesG;
    delete [] extractedValuesG;
    delete [] extractedIndicesS;
    delete [] extractedValuesS;
    return diag;
  }