1: #include <petscconf.h>
  2:  #include <../src/mat/impls/aij/mpi/mpiaij.h>
  3:  #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>

  5: PetscErrorCode  MatMPIAIJSetPreallocation_MPIAIJViennaCL(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[])
  6: {
  7:   Mat_MPIAIJ     *b = (Mat_MPIAIJ*)B->data;

 11:   PetscLayoutSetUp(B->rmap);
 12:   PetscLayoutSetUp(B->cmap);
 13:   if (!B->preallocated) {
 14:     /* Explicitly create the two MATSEQAIJVIENNACL matrices. */
 15:     MatCreate(PETSC_COMM_SELF,&b->A);
 16:     MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);
 17:     MatSetType(b->A,MATSEQAIJVIENNACL);
 18:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);
 19:     MatCreate(PETSC_COMM_SELF,&b->B);
 20:     MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);
 21:     MatSetType(b->B,MATSEQAIJVIENNACL);
 22:     PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);
 23:   }
 24:   MatSeqAIJSetPreallocation(b->A,d_nz,d_nnz);
 25:   MatSeqAIJSetPreallocation(b->B,o_nz,o_nnz);
 26:   B->preallocated = PETSC_TRUE;
 27:   return(0);
 28: }

 30: PetscErrorCode  MatCreateVecs_MPIAIJViennaCL(Mat mat,Vec *right,Vec *left)
 31: {
 33:   PetscInt       rbs,cbs;

 36:   MatGetBlockSizes(mat,&rbs,&cbs);
 37:   if (right) {
 38:     VecCreate(PetscObjectComm((PetscObject)mat),right);
 39:     VecSetSizes(*right,mat->cmap->n,PETSC_DETERMINE);
 40:     VecSetBlockSize(*right,cbs);
 41:     VecSetType(*right,VECVIENNACL);
 42:     VecSetLayout(*right,mat->cmap);
 43:   }
 44:   if (left) {
 45:     VecCreate(PetscObjectComm((PetscObject)mat),left);
 46:     VecSetSizes(*left,mat->rmap->n,PETSC_DETERMINE);
 47:     VecSetBlockSize(*left,rbs);
 48:     VecSetType(*left,VECVIENNACL);
 49:     VecSetLayout(*left,mat->rmap);
 50:   }
 51:   return(0);
 52: }

 54: PetscErrorCode MatAssemblyEnd_MPIAIJViennaCL(Mat A,MatAssemblyType mode)
 55: {
 56:   Mat_MPIAIJ     *b = (Mat_MPIAIJ*)A->data;
 58:   PetscBool      v;

 61:   MatAssemblyEnd_MPIAIJ(A,mode);
 62:   PetscObjectTypeCompare((PetscObject)b->lvec,VECSEQVIENNACL,&v);
 63:   if (!v) {
 64:     PetscInt m;
 65:     VecGetSize(b->lvec,&m);
 66:     VecDestroy(&b->lvec);
 67:     VecCreateSeqViennaCL(PETSC_COMM_SELF,m,&b->lvec);
 68:   }
 69:   return(0);
 70: }

 72: PetscErrorCode MatDestroy_MPIAIJViennaCL(Mat A)
 73: {

 77:   MatDestroy_MPIAIJ(A);
 78:   return(0);
 79: }

 81: PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJViennaCL(Mat A)
 82: {

 86:   MatCreate_MPIAIJ(A);
 87:   PetscObjectComposeFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",MatMPIAIJSetPreallocation_MPIAIJViennaCL);
 88:   A->ops->getvecs      = MatCreateVecs_MPIAIJViennaCL;
 89:   A->ops->assemblyend  = MatAssemblyEnd_MPIAIJViennaCL;
 90:   PetscObjectChangeTypeName((PetscObject)A,MATMPIAIJVIENNACL);
 91:   return(0);
 92: }


 95: /*@
 96:    MatCreateAIJViennaCL - Creates a sparse matrix in AIJ (compressed row) format
 97:    (the default parallel PETSc format).  This matrix will ultimately be pushed down
 98:    to GPUs and use the ViennaCL library for calculations. For good matrix
 99:    assembly performance the user should preallocate the matrix storage by setting
100:    the parameter nz (or the array nnz).  By setting these parameters accurately,
101:    performance during matrix assembly can be increased substantially.


104:    Collective on MPI_Comm

106:    Input Parameters:
107: +  comm - MPI communicator, set to PETSC_COMM_SELF
108: .  m - number of rows
109: .  n - number of columns
110: .  nz - number of nonzeros per row (same for all rows)
111: -  nnz - array containing the number of nonzeros in the various rows
112:          (possibly different for each row) or NULL

114:    Output Parameter:
115: .  A - the matrix

117:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
118:    MatXXXXSetPreallocation() paradigm instead of this routine directly.
119:    [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]

121:    Notes:
122:    If nnz is given then nz is ignored

124:    The AIJ format (also called the Yale sparse matrix format or
125:    compressed row storage), is fully compatible with standard Fortran 77
126:    storage.  That is, the stored row and column indices can begin at
127:    either one (as in Fortran) or zero.  See the users' manual for details.

129:    Specify the preallocated storage with either nz or nnz (not both).
130:    Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
131:    allocation.  For large problems you MUST preallocate memory or you
132:    will get TERRIBLE performance, see the users' manual chapter on matrices.

134:    Level: intermediate

136: .seealso: MatCreate(), MatCreateAIJ(), MatCreateAIJCUSPARSE(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatCreateAIJ(), MATMPIAIJVIENNACL, MATAIJVIENNACL
137: @*/
138: PetscErrorCode  MatCreateAIJViennaCL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[],Mat *A)
139: {
141:   PetscMPIInt    size;

144:   MatCreate(comm,A);
145:   MatSetSizes(*A,m,n,M,N);
146:   MPI_Comm_size(comm,&size);
147:   if (size > 1) {
148:     MatSetType(*A,MATMPIAIJVIENNACL);
149:     MatMPIAIJSetPreallocation(*A,d_nz,d_nnz,o_nz,o_nnz);
150:   } else {
151:     MatSetType(*A,MATSEQAIJVIENNACL);
152:     MatSeqAIJSetPreallocation(*A,d_nz,d_nnz);
153:   }
154:   return(0);
155: }

157: /*MC
158:    MATAIJVIENNACL - MATMPIAIJVIENNACL= "aijviennacl" = "mpiaijviennacl" - A matrix type to be used for sparse matrices.

160:    A matrix type (CSR format) whose data resides on GPUs.
161:    All matrix calculations are performed using the ViennaCL library.

163:    This matrix type is identical to MATSEQAIJVIENNACL when constructed with a single process communicator,
164:    and MATMPIAIJVIENNACL otherwise.  As a result, for single process communicators,
165:    MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported
166:    for communicators controlling multiple processes.  It is recommended that you call both of
167:    the above preallocation routines for simplicity.

169:    Options Database Keys:
170: +  -mat_type mpiaijviennacl - sets the matrix type to "mpiaijviennacl" during a call to MatSetFromOptions()

172:   Level: beginner

174:  .seealso: MatCreateAIJViennaCL(), MATSEQAIJVIENNACL, MatCreateSeqAIJVIENNACL()
175: M*/