mat creation

[Barry Smith]

On Mon, 6 Feb 2006, Roberto Gori wrote:

> Many thanks Barry,
>
> MatMPIAIJSetPreallocationCSR could be a good solution for me.
> It seems that to satisfy the MatMPIAIJSetPreallocationCSR  requirements I should convert the nz array from a relative to an absolute representation of the row offsets but ...
> I think this implies a not negligible drawback: my matrix is too large and offsets must be promoted to a long long int (64 bit) type.
> Moreover the  idxn array type also must be promoted to be compliant with a 64 bit PetscInt  type.  This requires the double amount of memory for nz and idxn.
>
> What do you think about the opportunity to provide another version of  MatMPIAIJSetPreallocationCSR  based on relative offsets to save memory?
> Absolute offsets could be managed internally  ...

  Your point is valid and noted but this would require writing another
whole set of the AIJ routines that use the relative offsets. Maybe not so
hard but a lot of code duplication.

>
> Finally I've look into the body of  MatMPIAIJSetPreallocationCSR_MPIAIJ function ad I found this snippet of code
>
> for (i=0; i<m; i++) {
>     ii   = i + rstart;
>     nnz  = I[i+1]- I[i];
>     ierr = MatSetValues_MPIAIJ(B,1,&ii,nnz,J+I[i],values+(v ? I[i] : 0),INSERT_VALUES);CHKERRQ(ierr);
>   }
>
> There's no way to avoid the copy of the values? This doubles one more time the memory for the values array and spend time into the loop.
> Morevover nnz corresponds exactly to my nz[i] ...

  The thing is PETSc does not just store the MPIAIJ matrices in a simple
CSR format (the matrix is actually stored in two parts, see the MATMPIAIJ
manual page). So the copy is needed.

  What we do. Is first count the number of nonzeros per row then
preallocate the memory then compute the values in the matrix, see
MatMPIAIJSetPreallocatio(). This may
seem expensive but really is not, the counting is fast. You build the CSR
matrix for your sparse matrix; using our approach is no harder and
prevents the need for a copy. (In fact the reason PETSc has the
Preallocation() and MatSetValues() stuff is so that users do not have
to construct the CSR format directly.

   Barry



 >
> ----- Original Message -----
>   From: Barry Smith
>   To: Roberto Gori
>   Cc: petsc-users at mcs.anl.gov
>   Sent: Monday, February 06, 2006 4:15 PM
>   Subject: Re: mat creation
>
>
>
>
>   On Mon, 6 Feb 2006, Roberto Gori wrote:
>
>   > Hi, guys,
>   >
>   > I'm trying to solve a linear system with M equations and N unknowns (M >> N) using the LSQR method.
>   >
>   > I have a C parallel sparse matrix MxN with NZ nonzero elements stored in this way:
>   >
>   > double values[NZ]; // the nonzero values
>   > int idxn[NZ]; // the column indices
>   > int nz[M]; // the number of nonzero values for each row
>
>   This needs to be instead indices to start of each row in idxn and values
>   instead of the length of each row.
>   >
>   > what is the best way to create a PETSc matrix?
>
>   MatMPIAIJSetPreallocationCSR()
>
>
>   >
>   > Do exists a statement like VecCreateMPIWithArray for already distributed vectors that could do all the work in one shot without copy?
>
>   Yes
>
>   >  Just another detail. Initially M is unknown.  All my structures are
>   > allocated dynamically because the number of rows M (and the number
>   > of nonzero NZ also) is determined just before to solve the system,
>   > so I think it's no feasible for me to call MatCreateMPIAIJ at the
>   > begin of my code and to perform the MatSetValues loop.
>
>     There is nothing saying that you must "call MatCreateMPIAIJ at the
>   beginning of your code", you can call it anytime that you have the
>   needed information.
>
>      Barry
>
>   >
>   >     Roberto
>   >
>   > ================================================================
>   > | Roberto Gori                                            CINECA
>   > | High Performance System Group         Via Magnanelli 6/3 40033
>   > | e-mail: r.gori at cineca.it                          Casalecchio di Reno
>   > | Tel: 051/6171522                                    Bologna - ITALY
>   > | Fax: 051/6132198                                   http://www.cineca.it
>   > ================================================================
>   >
>   >
>
>