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