dolfin team mailing list archive

[jesperc <jesperc@xxxxxx>]

I found the problem with the segmentation fault. The vals array should 
contain all elements of the matrix which seems really inefficient for an 
identity matrix:

int n = 50000;
real *vals = new real[n*n];
int *cols = new int[n];
int *rows = new int[n];
for (int i=0; i<n; ++i)
 {
 for (int j=0; j<n; ++j)
   {
   if (i=j)
     vals[j+i*n] = 1.0;
   else
     vals[j+i*n] = 0.0;
   }
 rows[i] = i;
 cols[i] = i;
 }
A->set(vals, cols, n, rows, n);
A->apply();

Another way is to use set() for each row

int n = 50000;
real *vals = new real[n];
int *cols = new int[n];
int *rows = new int[1];
for (int i=0; i<n; ++i)
 {
 vals[i] = 0.0;
 cols[i] = i;
 }
for (int i=0; i<n; ++i)
 {
 rows[0] = i;
 vals[i] = 0.0;
 A->set(vals, cols, n, rows, 1);
 vals[i] = 0.0;
 }
A->apply();

Both ways seem to be really slow though.

/Jesper

Garth N. Wells skrev:
> jesperc wrote:
> > Hi again,
> >
> > I realize that assigning values blockwise is more efficient but if I 
> > try the following I get a segmentation fault (it works for one single 
> > row/column at a time though):
> >
> > Matrix A;
> > A.init(n,n,1);
> > real *vals = new real[n];
> > int *cols = new int[n];
> > int *rows = new int[n];
> > for (int i=0; i<n; ++i)
> >   {
> >    vals[i] = 1.0;
> >    cols[i] = i;
> >    rows[i] = i;
> >   }
> > A->set(vals, cols, n, rows, n);
> > A->apply();
>
>
> I don't know what's going on. You could try
>
>   A.init(n, n, 2);
>
> If you can test against the development version of DOLFIN and still 
> have a problem, I can take a look at it.
>
> Garth
>
> > Regards, Jesper
> >
> > Garth N. Wells skrev:
> > > jesperc wrote:
> > > > Hi,
> > > >
> > > > Does anyone have an example of using the ident() function in 
> > > > PETScMatrix. I cannot get it to work 
> > > Matrix::ident(. . . ) is not meant to do this, so it is unlikely to 
> > > work.
> > >
> > > and using a for loop with set() to
> > > > create an identity matrix is painfully slow...
> > > >
> > > > This is how I do it now:
> > > >
> > > > Matrix A;
> > > > A.init(n,n,1);
> > > > for (int i=0; i<n; ++i)
> > > >    A.set(i,i,1);
> > > You shouldn't set Matrix values element-wise 
> > > (Matrix::set(int,int,real) has been removed from the Matrix 
> > > interface to stop users doing this). If you really want to access 
> > > matrix entries element-wise, use a uBlasMatrix, otherwise set values 
> > > in large blocks using
> > >
> > >   Matrix::set(const real* block, uint m, const uint* rows, uint n, 
> > > const uint* cols)
> > >
> > > Garth
> > >
> > > > and it kills me/my computer.
> > > >
> > > > Kind regards, Jesper
> > > > _______________________________________________
> > > > DOLFIN-dev mailing list
> > > > DOLFIN-dev@xxxxxxxxxx
> > > > http://www.fenics.org/mailman/listinfo/dolfin-dev
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/mailman/listinfo/dolfin-dev