dolfin team mailing list archive

[Matthew Knepley <knepley@xxxxxxxxxxx>]

Anders Logg <logg@xxxxxxxxx> writes:

> I went with the simple option of using MatGetRow(), see below. 
> (Don't need the extra overhead of the parallel version right now.)

  Okay.

> Doing all scalar products by hand one by one instead of calling
> MatMult() on all rows at once is about a factor 3.5 slower (and it
> seems to get better when I refine the mesh).

  Thats better than I would have guesses.

> By the way, why do you have int and not unsigned int for size types,
> column arrays etc? I have to do casts all over the place. (We're also
> using STL which uses unsigned int and it propagates throughout the
> code.)

  We use negative numbers to indicate different conditions all the time.

     Thanks,

        Matt

> /Anders
>
> real Matrix::mult(const Vector& x, uint row) const
> {
>   // FIXME: Temporary fix (assumes uniprocessor case)
>
>   int ncols = 0;
>   const int* cols = 0;
>   const double* Avals = 0;
>   double* xvals = 0;
>   MatGetRow(A, static_cast<int>(row), &ncols, &cols, &Avals);
>   VecGetArray(x.x, &xvals);
>
>   real sum = 0.0;
>   for (int i = 0; i < ncols; i++)
>     sum += Avals[i] * xvals[cols[i]];
>
>   MatRestoreRow(A, static_cast<int>(row), &ncols, &cols, &Avals);
>   VecRestoreArray(x.x, &xvals);
>
>   return sum;
> }
>
> On Tue, Jun 07, 2005 at 08:37:37AM -0500, Barry Smith wrote:
>> 
>>    Not again, Barry
>> 
>> On Tue, 7 Jun 2005, Barry Smith wrote:
>> 
>> > 
>> >   Missed one thing,
>> > 
>> > >   I disagree. For instance, we have the MatGetLocalMat() function. Further, I
>> > > think that this can be generalized. Most parallel formats will have something
>> >   ^^^^^^
>> > 
>> >    I prefer never to generalize until I KNOW it can be generalized.
>> 
>>    and that the generalization actually buys you something.
>> 
>>   Barry
>> 
>> > 
>> >    Barry
>> > 
>> > On Tue, 7 Jun 2005, Matthew Knepley wrote:
>> > 
>> > > Barry Smith <petsc-maint@xxxxxxxxxxx> writes:
>> > > 
>> > > >   Matt,
>> > > >
>> > > >   If you are coding it specifically for MPIAIJ (as you have done) why introduce
>> > > > the really ugly MatGetCommunicationStructs() (which is very specific to certain
>> > > > matrix formats anyways) just access the scatter info directly and bag that 
>> > > > function. 
>> > > 
>> > >   I disagree. For instance, we have the MatGetLocalMat() function. Further, I
>> > > think that this can be generalized. Most parallel formats will have something
>> > > like this.
>> > > 
>> > >         Matt
>> > > 
>> > > >    Barry
>> > > >
>> > > >
>> > > > On Tue, 7 Jun 2005, Matthew Knepley wrote:
>> > > >
>> > > >> Anders Logg <logg@xxxxxxxxx> writes:
>> > > >> 
>> > > >>        Alright, I have made the updates to PETSc. I did not add MatMultRow()
>> > > >> since there are so many efficiency implications, that it is best left to the
>> > > >> user. Here is my prototype code:
>> > > >> 
>> > > >> #undef __FUNCT__  
>> > > >> #define __FUNCT__ "MatMultRowMPIAIJ"
>> > > >> PetscErrorCode MatMultRowMPIAIJ(Mat A,Vec xx,int row,double *yy)
>> > > >> {
>> > > >>   Vec            lvec;
>> > > >> #if defined (PETSC_USE_CTABLE)
>> > > >>   PetscTable     colmap;
>> > > >> #else
>> > > >>   PetscInt      *colmap;
>> > > >> #endif
>> > > >>   VecScatter     mulScatter;
>> > > >>   PetscInt       ncols, lcol, start, end;
>> > > >>   PetscInt      *cols;
>> > > >>   PetscScalar   *vals, *array, *larray;
>> > > >>   PetscInt       nt, i;
>> > > >>   PetscErrorCode ierr;
>> > > >> 
>> > > >>   PetscFunctionBegin;
>> > > >>   ierr = VecGetLocalSize(xx,&nt);CHKERRQ(ierr);
>> > > >>   if (nt != A->n) {
>> > > >>     SETERRQ2(PETSC_ERR_ARG_SIZ,"Incompatible partition of A (%D) and xx (%D)",A->n,nt);
>> > > >>   }
>> > > >>   ierr = MatGetCommunicationStructs(A, &lvec, &colmap, &multScatter);CHKERRQ(ierr);
>> > > >>   ierr = VecScatterBegin(xx,lvec,INSERT_VALUES,SCATTER_FORWARD,multScatter);CHKERRQ(ierr);
>> > > >>   ierr = VecScatterEnd(xx,lvec,INSERT_VALUES,SCATTER_FORWARD,multScatter);CHKERRQ(ierr);
>> > > >>   ierr = VecGetArray(xx, &array);CHKERRQ(ierr);
>> > > >>   ierr = VecGetArray(lvec, &larray);CHKERRQ(ierr);
>> > > >>   ierr = MatGetOwnershipRange(A, &start, &end);CHKERRQ(ierr);
>> > > >>   ierr = MatGetRow(A, row, &ncols, &cols, &vals);CHKERRQ(ierr);
>> > > >>   *yy = 0.0;
>> > > >>   for(i = 0; i < ncols; i++) {
>> > > >>     if ((col[i] >= start) && (col[i] < end)) {
>> > > >>       *yy += vals[i]*array[col[i]-start];
>> > > >>     } else {
>> > > >> #if defined (PETSC_USE_CTABLE)
>> > > >>       ierr = PetscTableFind(colmap,col[i],&lcol);CHKERRQ(ierr);
>> > > >> #else
>> > > >>       lcol = colmap[col[i]];
>> > > >> #endif
>> > > >>       *yy += vals[i]*larray[lcol];
>> > > >>     }
>> > > >>   }
>> > > >>   ierr = MatRestoreRow(A, row, &ncols, &cols, &vals);CHKERRQ(ierr);
>> > > >>   ierr = VecRestoreArray(xx, &array);CHKERRQ(ierr);
>> > > >>   ierr = VecRestoreArray(lvec, &larray);CHKERRQ(ierr);
>> > > >>   PetscFunctionReturn(0);
>> > > >> }
>> > > >> 
>> > > >> You just have to pull and rebuild in src/mat.
>> > > >> 
>> > > >>         Thanks,
>> > > >> 
>> > > >>                 Matt
>> > > >> 
>> > > >> > Thanks. Please take a look if you have time since it is something I
>> > > >> > really need. (Doing the whole multiply each time I need the value for
>> > > >> > one component would cost too much, and I need the different components
>> > > >> > at different times.)
>> > > >> >
>> > > >> > In the meantime, I will make a copy of the values and do the
>> > > >> > multiplication by hand.
>> > > >> >
>> > > >> > /Anders
>> > > >> >
>> > > >> > On Mon, Jun 06, 2005 at 10:06:23AM -0500, Matthew Knepley wrote:
>> > > >> >> Anders Logg <logg@xxxxxxxxx> writes:
>> > > >> >> 
>> > > >> >> > Matt, how do I compute the scalar product of a row of a matrix and a
>> > > >> >> > given vector in PETSc?
>> > > >> >> >
>> > > >> >> > I need something like
>> > > >> >> >
>> > > >> >> >     MatMult(Mat A, Vec x, int i, double* yi)
>> > > >> >> >
>> > > >> >> > where yi = (A*x)_i for my multi-adaptive solver (need to evaluate each
>> > > >> >> > component of the right-hand side f individually).
>> > > >> >> 
>> > > >> >>   I don't think we have anything like that. Of course, you can just pick that
>> > > >> >> component out of the result vector, but I am guessing you do not want to do
>> > > >> >> the whole multiply. I will look at it, but my guess is that it would involve
>> > > >> >> a good bit of rewriting.
>> > > >> >> 
>> > > >> >>   Thanks,
>> > > >> >> 
>> > > >> >>         Matt
>> > > >> >
>> > > >> >
>> > > >> >
>> > > >> >
>> > > >> 
>> > > >> 
>> > > >
>> > > >
>> > > >
>> > > 
>> > > 
>> > 
>> 
>
> -- 
> Anders Logg
> Research Assistant Professor
> Toyota Technological Institute at Chicago
> http://www.tti-c.org/logg/
>
>
>

-- 
"Failure has a thousand explanations. Success doesn't need one" -- Sir Alec Guiness