xref: /petsc/src/mat/impls/aij/mpi/mpiviennacl/mpiaijviennacl.cxx (revision 37eeb8152ec6a2cf24186d3591c2c5de5dfd8fa5)
1 #include <petscconf.h>
2 #include <../src/mat/impls/aij/mpi/mpiaij.h>   /*I "petscmat.h" I*/
3 #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
4 
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;
8   PetscErrorCode ierr;
9 
10   PetscFunctionBegin;
11   ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr);
12   ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr);
13   if (!B->preallocated) {
14     /* Explicitly create the two MATSEQAIJVIENNACL matrices. */
15     ierr = MatCreate(PETSC_COMM_SELF,&b->A);CHKERRQ(ierr);
16     ierr = MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n);CHKERRQ(ierr);
17     ierr = MatSetType(b->A,MATSEQAIJVIENNACL);CHKERRQ(ierr);
18     ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)b->A);CHKERRQ(ierr);
19     ierr = MatCreate(PETSC_COMM_SELF,&b->B);CHKERRQ(ierr);
20     ierr = MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N);CHKERRQ(ierr);
21     ierr = MatSetType(b->B,MATSEQAIJVIENNACL);CHKERRQ(ierr);
22     ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)b->B);CHKERRQ(ierr);
23   }
24   ierr = MatSeqAIJSetPreallocation(b->A,d_nz,d_nnz);CHKERRQ(ierr);
25   ierr = MatSeqAIJSetPreallocation(b->B,o_nz,o_nnz);CHKERRQ(ierr);
26   B->preallocated = PETSC_TRUE;
27   PetscFunctionReturn(0);
28 }
29 
30 PetscErrorCode MatAssemblyEnd_MPIAIJViennaCL(Mat A,MatAssemblyType mode)
31 {
32   Mat_MPIAIJ     *b = (Mat_MPIAIJ*)A->data;
33   PetscErrorCode ierr;
34   PetscBool      v;
35 
36   PetscFunctionBegin;
37   ierr = MatAssemblyEnd_MPIAIJ(A,mode);CHKERRQ(ierr);
38   ierr = PetscObjectTypeCompare((PetscObject)b->lvec,VECSEQVIENNACL,&v);CHKERRQ(ierr);
39   if (!v) {
40     PetscInt m;
41     ierr = VecGetSize(b->lvec,&m);CHKERRQ(ierr);
42     ierr = VecDestroy(&b->lvec);CHKERRQ(ierr);
43     ierr = VecCreateSeqViennaCL(PETSC_COMM_SELF,m,&b->lvec);CHKERRQ(ierr);
44   }
45   PetscFunctionReturn(0);
46 }
47 
48 PetscErrorCode MatDestroy_MPIAIJViennaCL(Mat A)
49 {
50   PetscErrorCode ierr;
51 
52   PetscFunctionBegin;
53   ierr = MatDestroy_MPIAIJ(A);CHKERRQ(ierr);
54   PetscFunctionReturn(0);
55 }
56 
57 PETSC_EXTERN PetscErrorCode MatCreate_MPIAIJViennaCL(Mat A)
58 {
59   PetscErrorCode ierr;
60 
61   PetscFunctionBegin;
62   ierr = MatCreate_MPIAIJ(A);CHKERRQ(ierr);
63   ierr = PetscFree(A->defaultvectype);CHKERRQ(ierr);
64   ierr = PetscStrallocpy(VECVIENNACL,&A->defaultvectype);CHKERRQ(ierr);
65   ierr = PetscObjectComposeFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",MatMPIAIJSetPreallocation_MPIAIJViennaCL);CHKERRQ(ierr);
66   A->ops->assemblyend = MatAssemblyEnd_MPIAIJViennaCL;
67   ierr = PetscObjectChangeTypeName((PetscObject)A,MATMPIAIJVIENNACL);CHKERRQ(ierr);
68   PetscFunctionReturn(0);
69 }
70 
71 
72 /*@
73    MatCreateAIJViennaCL - Creates a sparse matrix in AIJ (compressed row) format
74    (the default parallel PETSc format).  This matrix will ultimately be pushed down
75    to GPUs and use the ViennaCL library for calculations. For good matrix
76    assembly performance the user should preallocate the matrix storage by setting
77    the parameter nz (or the array nnz).  By setting these parameters accurately,
78    performance during matrix assembly can be increased substantially.
79 
80 
81    Collective
82 
83    Input Parameters:
84 +  comm - MPI communicator, set to PETSC_COMM_SELF
85 .  m - number of rows
86 .  n - number of columns
87 .  nz - number of nonzeros per row (same for all rows)
88 -  nnz - array containing the number of nonzeros in the various rows
89          (possibly different for each row) or NULL
90 
91    Output Parameter:
92 .  A - the matrix
93 
94    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
95    MatXXXXSetPreallocation() paradigm instead of this routine directly.
96    [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
97 
98    Notes:
99    If nnz is given then nz is ignored
100 
101    The AIJ format (also called the Yale sparse matrix format or
102    compressed row storage), is fully compatible with standard Fortran 77
103    storage.  That is, the stored row and column indices can begin at
104    either one (as in Fortran) or zero.  See the users' manual for details.
105 
106    Specify the preallocated storage with either nz or nnz (not both).
107    Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
108    allocation.  For large problems you MUST preallocate memory or you
109    will get TERRIBLE performance, see the users' manual chapter on matrices.
110 
111    Level: intermediate
112 
113 .seealso: MatCreate(), MatCreateAIJ(), MatCreateAIJCUSPARSE(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatCreateAIJ(), MATMPIAIJVIENNACL, MATAIJVIENNACL
114 @*/
115 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)
116 {
117   PetscErrorCode ierr;
118   PetscMPIInt    size;
119 
120   PetscFunctionBegin;
121   ierr = MatCreate(comm,A);CHKERRQ(ierr);
122   ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr);
123   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
124   if (size > 1) {
125     ierr = MatSetType(*A,MATMPIAIJVIENNACL);CHKERRQ(ierr);
126     ierr = MatMPIAIJSetPreallocation(*A,d_nz,d_nnz,o_nz,o_nnz);CHKERRQ(ierr);
127   } else {
128     ierr = MatSetType(*A,MATSEQAIJVIENNACL);CHKERRQ(ierr);
129     ierr = MatSeqAIJSetPreallocation(*A,d_nz,d_nnz);CHKERRQ(ierr);
130   }
131   PetscFunctionReturn(0);
132 }
133 
134 /*MC
135    MATAIJVIENNACL - MATMPIAIJVIENNACL= "aijviennacl" = "mpiaijviennacl" - A matrix type to be used for sparse matrices.
136 
137    A matrix type (CSR format) whose data resides on GPUs.
138    All matrix calculations are performed using the ViennaCL library.
139 
140    This matrix type is identical to MATSEQAIJVIENNACL when constructed with a single process communicator,
141    and MATMPIAIJVIENNACL otherwise.  As a result, for single process communicators,
142    MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported
143    for communicators controlling multiple processes.  It is recommended that you call both of
144    the above preallocation routines for simplicity.
145 
146    Options Database Keys:
147 +  -mat_type mpiaijviennacl - sets the matrix type to "mpiaijviennacl" during a call to MatSetFromOptions()
148 
149   Level: beginner
150 
151  .seealso: MatCreateAIJViennaCL(), MATSEQAIJVIENNACL, MatCreateSeqAIJVIENNACL()
152 M*/
153