xref: /petsc/src/mat/impls/aij/seq/seqviennacl/aijviennacl.cxx (revision 5e71baeff2f3138f93cd4f5927dfd596eb8325cc)
1 
2 
3 /*
4     Defines the basic matrix operations for the AIJ (compressed row)
5   matrix storage format.
6 */
7 
8 #include <petscconf.h>
9 #include <../src/mat/impls/aij/seq/aij.h>          /*I "petscmat.h" I*/
10 #include <petscbt.h>
11 #include <../src/vec/vec/impls/dvecimpl.h>
12 #include <petsc/private/vecimpl.h>
13 
14 #include <../src/mat/impls/aij/seq/seqviennacl/viennaclmatimpl.h>
15 
16 
17 #include <algorithm>
18 #include <vector>
19 #include <string>
20 
21 #include "viennacl/linalg/prod.hpp"
22 
23 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJViennaCL(Mat A, MatType type, MatReuse reuse, Mat *newmat);
24 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*);
25 
26 
27 PetscErrorCode MatViennaCLCopyToGPU(Mat A)
28 {
29 
30   Mat_SeqAIJViennaCL *viennaclstruct = (Mat_SeqAIJViennaCL*)A->spptr;
31   Mat_SeqAIJ         *a              = (Mat_SeqAIJ*)A->data;
32   PetscErrorCode     ierr;
33 
34 
35   PetscFunctionBegin;
36   if (A->rmap->n > 0 && A->cmap->n > 0) { //some OpenCL SDKs have issues with buffers of size 0
37     if (A->valid_GPU_matrix == PETSC_VIENNACL_UNALLOCATED || A->valid_GPU_matrix == PETSC_VIENNACL_CPU) {
38       ierr = PetscLogEventBegin(MAT_ViennaCLCopyToGPU,A,0,0,0);CHKERRQ(ierr);
39 
40       try {
41         if (a->compressedrow.use) {
42           if (!viennaclstruct->compressed_mat) viennaclstruct->compressed_mat = new ViennaCLCompressedAIJMatrix();
43 
44           // Since PetscInt is different from cl_uint, we have to convert:
45           viennacl::backend::mem_handle dummy;
46 
47           viennacl::backend::typesafe_host_array<unsigned int> row_buffer; row_buffer.raw_resize(dummy, a->compressedrow.nrows+1);
48           for (PetscInt i=0; i<=a->compressedrow.nrows; ++i)
49             row_buffer.set(i, (a->compressedrow.i)[i]);
50 
51           viennacl::backend::typesafe_host_array<unsigned int> row_indices; row_indices.raw_resize(dummy, a->compressedrow.nrows);
52           for (PetscInt i=0; i<a->compressedrow.nrows; ++i)
53             row_indices.set(i, (a->compressedrow.rindex)[i]);
54 
55           viennacl::backend::typesafe_host_array<unsigned int> col_buffer; col_buffer.raw_resize(dummy, a->nz);
56           for (PetscInt i=0; i<a->nz; ++i)
57             col_buffer.set(i, (a->j)[i]);
58 
59           viennaclstruct->compressed_mat->set(row_buffer.get(), row_indices.get(), col_buffer.get(), a->a, A->rmap->n, A->cmap->n, a->compressedrow.nrows, a->nz);
60         } else {
61           if (!viennaclstruct->mat) viennaclstruct->mat = new ViennaCLAIJMatrix();
62 
63           // Since PetscInt is in general different from cl_uint, we have to convert:
64           viennacl::backend::mem_handle dummy;
65 
66           viennacl::backend::typesafe_host_array<unsigned int> row_buffer; row_buffer.raw_resize(dummy, A->rmap->n+1);
67           for (PetscInt i=0; i<=A->rmap->n; ++i)
68             row_buffer.set(i, (a->i)[i]);
69 
70           viennacl::backend::typesafe_host_array<unsigned int> col_buffer; col_buffer.raw_resize(dummy, a->nz);
71           for (PetscInt i=0; i<a->nz; ++i)
72             col_buffer.set(i, (a->j)[i]);
73 
74           viennaclstruct->mat->set(row_buffer.get(), col_buffer.get(), a->a, A->rmap->n, A->cmap->n, a->nz);
75         }
76         ViennaCLWaitForGPU();
77       } catch(std::exception const & ex) {
78         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex.what());
79       }
80 
81       // Create temporary vector for v += A*x:
82       if (viennaclstruct->tempvec) {
83         if (viennaclstruct->tempvec->size() != static_cast<std::size_t>(A->rmap->n)) {
84           delete (ViennaCLVector*)viennaclstruct->tempvec;
85           viennaclstruct->tempvec = new ViennaCLVector(A->rmap->n);
86         } else {
87           viennaclstruct->tempvec->clear();
88         }
89       } else {
90         viennaclstruct->tempvec = new ViennaCLVector(A->rmap->n);
91       }
92 
93       A->valid_GPU_matrix = PETSC_VIENNACL_BOTH;
94 
95       ierr = PetscLogEventEnd(MAT_ViennaCLCopyToGPU,A,0,0,0);CHKERRQ(ierr);
96     }
97   }
98   PetscFunctionReturn(0);
99 }
100 
101 PetscErrorCode MatViennaCLCopyFromGPU(Mat A, const ViennaCLAIJMatrix *Agpu)
102 {
103   Mat_SeqAIJ         *a              = (Mat_SeqAIJ*)A->data;
104   PetscInt           m               = A->rmap->n;
105   PetscErrorCode     ierr;
106 
107 
108   PetscFunctionBegin;
109   if (A->valid_GPU_matrix == PETSC_VIENNACL_UNALLOCATED) {
110     try {
111       if (a->compressedrow.use) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "ViennaCL: Cannot handle row compression for GPU matrices");
112       else {
113 
114         if ((PetscInt)Agpu->size1() != m) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "GPU matrix has %d rows, should be %d", Agpu->size1(), m);
115         a->nz           = Agpu->nnz();
116         a->maxnz        = a->nz; /* Since we allocate exactly the right amount */
117         A->preallocated = PETSC_TRUE;
118         if (a->singlemalloc) {
119           if (a->a) {ierr = PetscFree3(a->a,a->j,a->i);CHKERRQ(ierr);}
120         } else {
121           if (a->i) {ierr = PetscFree(a->i);CHKERRQ(ierr);}
122           if (a->j) {ierr = PetscFree(a->j);CHKERRQ(ierr);}
123           if (a->a) {ierr = PetscFree(a->a);CHKERRQ(ierr);}
124         }
125         ierr = PetscMalloc3(a->nz,&a->a,a->nz,&a->j,m+1,&a->i);CHKERRQ(ierr);
126         ierr = PetscLogObjectMemory((PetscObject)A, a->nz*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr);
127 
128         a->singlemalloc = PETSC_TRUE;
129 
130         /* Setup row lengths */
131         if (a->imax) {ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr);}
132         ierr = PetscMalloc2(m,&a->imax,m,&a->ilen);CHKERRQ(ierr);
133         ierr = PetscLogObjectMemory((PetscObject)A, 2*m*sizeof(PetscInt));CHKERRQ(ierr);
134 
135         /* Copy data back from GPU */
136         viennacl::backend::typesafe_host_array<unsigned int> row_buffer; row_buffer.raw_resize(Agpu->handle1(), Agpu->size1() + 1);
137 
138         // copy row array
139         viennacl::backend::memory_read(Agpu->handle1(), 0, row_buffer.raw_size(), row_buffer.get());
140         (a->i)[0] = row_buffer[0];
141         for (PetscInt i = 0; i < (PetscInt)Agpu->size1(); ++i) {
142           (a->i)[i+1] = row_buffer[i+1];
143           a->imax[i]  = a->ilen[i] = a->i[i+1] - a->i[i];  //Set imax[] and ilen[] arrays at the same time as i[] for better cache reuse
144         }
145 
146         // copy column indices
147         viennacl::backend::typesafe_host_array<unsigned int> col_buffer; col_buffer.raw_resize(Agpu->handle2(), Agpu->nnz());
148         viennacl::backend::memory_read(Agpu->handle2(), 0, col_buffer.raw_size(), col_buffer.get());
149         for (PetscInt i=0; i < (PetscInt)Agpu->nnz(); ++i)
150           (a->j)[i] = col_buffer[i];
151 
152         // copy nonzero entries directly to destination (no conversion required)
153         viennacl::backend::memory_read(Agpu->handle(), 0, sizeof(PetscScalar)*Agpu->nnz(), a->a);
154 
155         ViennaCLWaitForGPU();
156         /* TODO: Once a->diag is moved out of MatAssemblyEnd(), invalidate it here. */
157       }
158     } catch(std::exception const & ex) {
159       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex.what());
160     }
161 
162     /* This assembly prevents resetting the flag to PETSC_VIENNACL_CPU and recopying */
163     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
164     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
165 
166     A->valid_GPU_matrix = PETSC_VIENNACL_BOTH;
167   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "ViennaCL error: Only valid for unallocated GPU matrices");
168   PetscFunctionReturn(0);
169 }
170 
171 PetscErrorCode MatCreateVecs_SeqAIJViennaCL(Mat mat, Vec *right, Vec *left)
172 {
173   PetscErrorCode ierr;
174   PetscInt rbs,cbs;
175 
176   PetscFunctionBegin;
177   ierr = MatGetBlockSizes(mat,&rbs,&cbs);CHKERRQ(ierr);
178   if (right) {
179     ierr = VecCreate(PetscObjectComm((PetscObject)mat),right);CHKERRQ(ierr);
180     ierr = VecSetSizes(*right,mat->cmap->n,PETSC_DETERMINE);CHKERRQ(ierr);
181     ierr = VecSetBlockSize(*right,cbs);CHKERRQ(ierr);
182     ierr = VecSetType(*right,VECSEQVIENNACL);CHKERRQ(ierr);
183     ierr = PetscLayoutReference(mat->cmap,&(*right)->map);CHKERRQ(ierr);
184   }
185   if (left) {
186     ierr = VecCreate(PetscObjectComm((PetscObject)mat),left);CHKERRQ(ierr);
187     ierr = VecSetSizes(*left,mat->rmap->n,PETSC_DETERMINE);CHKERRQ(ierr);
188     ierr = VecSetBlockSize(*left,rbs);CHKERRQ(ierr);
189     ierr = VecSetType(*left,VECSEQVIENNACL);CHKERRQ(ierr);
190     ierr = PetscLayoutReference(mat->rmap,&(*left)->map);CHKERRQ(ierr);
191   }
192   PetscFunctionReturn(0);
193 }
194 
195 PetscErrorCode MatMult_SeqAIJViennaCL(Mat A,Vec xx,Vec yy)
196 {
197   Mat_SeqAIJ           *a = (Mat_SeqAIJ*)A->data;
198   PetscErrorCode       ierr;
199   Mat_SeqAIJViennaCL   *viennaclstruct = (Mat_SeqAIJViennaCL*)A->spptr;
200   const ViennaCLVector *xgpu=NULL;
201   ViennaCLVector       *ygpu=NULL;
202 
203   PetscFunctionBegin;
204   if (A->rmap->n > 0 && A->cmap->n > 0) {
205     ierr = VecViennaCLGetArrayRead(xx,&xgpu);CHKERRQ(ierr);
206     ierr = VecViennaCLGetArrayWrite(yy,&ygpu);CHKERRQ(ierr);
207     try {
208       *ygpu = viennacl::linalg::prod(*viennaclstruct->mat,*xgpu);
209       ViennaCLWaitForGPU();
210     } catch (std::exception const & ex) {
211       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex.what());
212     }
213     ierr = VecViennaCLRestoreArrayRead(xx,&xgpu);CHKERRQ(ierr);
214     ierr = VecViennaCLRestoreArrayWrite(yy,&ygpu);CHKERRQ(ierr);
215     ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr);
216   }
217   PetscFunctionReturn(0);
218 }
219 
220 
221 
222 PetscErrorCode MatMultAdd_SeqAIJViennaCL(Mat A,Vec xx,Vec yy,Vec zz)
223 {
224   Mat_SeqAIJ           *a = (Mat_SeqAIJ*)A->data;
225   PetscErrorCode       ierr;
226   Mat_SeqAIJViennaCL   *viennaclstruct = (Mat_SeqAIJViennaCL*)A->spptr;
227   const ViennaCLVector *xgpu=NULL,*ygpu=NULL;
228   ViennaCLVector       *zgpu=NULL;
229 
230   PetscFunctionBegin;
231   if (A->rmap->n > 0 && A->cmap->n > 0) {
232     try {
233       ierr = VecViennaCLGetArrayRead(xx,&xgpu);CHKERRQ(ierr);
234       ierr = VecViennaCLGetArrayRead(yy,&ygpu);CHKERRQ(ierr);
235       ierr = VecViennaCLGetArrayWrite(zz,&zgpu);CHKERRQ(ierr);
236 
237       if (a->compressedrow.use) {
238         ViennaCLVector temp = viennacl::linalg::prod(*viennaclstruct->compressed_mat, *xgpu);
239         *zgpu = *ygpu + temp;
240         ViennaCLWaitForGPU();
241       } else {
242         if (zz == xx || zz == yy) { //temporary required
243           ViennaCLVector temp = viennacl::linalg::prod(*viennaclstruct->mat, *xgpu);
244           *zgpu = *ygpu;
245           *zgpu += temp;
246           ViennaCLWaitForGPU();
247         } else {
248           *viennaclstruct->tempvec = viennacl::linalg::prod(*viennaclstruct->mat, *xgpu);
249           *zgpu = *ygpu + *viennaclstruct->tempvec;
250           ViennaCLWaitForGPU();
251         }
252       }
253 
254       ierr = VecViennaCLRestoreArrayRead(xx,&xgpu);CHKERRQ(ierr);
255       ierr = VecViennaCLRestoreArrayRead(yy,&ygpu);CHKERRQ(ierr);
256       ierr = VecViennaCLRestoreArrayWrite(zz,&zgpu);CHKERRQ(ierr);
257 
258     } catch(std::exception const & ex) {
259       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex.what());
260     }
261     ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr);
262   }
263   PetscFunctionReturn(0);
264 }
265 
266 PetscErrorCode MatAssemblyEnd_SeqAIJViennaCL(Mat A,MatAssemblyType mode)
267 {
268   PetscErrorCode ierr;
269 
270   PetscFunctionBegin;
271   ierr = MatAssemblyEnd_SeqAIJ(A,mode);CHKERRQ(ierr);
272   ierr = MatViennaCLCopyToGPU(A);CHKERRQ(ierr);
273   PetscFunctionReturn(0);
274 }
275 
276 /* --------------------------------------------------------------------------------*/
277 /*@
278    MatCreateSeqAIJViennaCL - Creates a sparse matrix in AIJ (compressed row) format
279    (the default parallel PETSc format).  This matrix will ultimately be pushed down
280    to GPUs and use the ViennaCL library for calculations. For good matrix
281    assembly performance the user should preallocate the matrix storage by setting
282    the parameter nz (or the array nnz).  By setting these parameters accurately,
283    performance during matrix assembly can be increased substantially.
284 
285 
286    Collective on MPI_Comm
287 
288    Input Parameters:
289 +  comm - MPI communicator, set to PETSC_COMM_SELF
290 .  m - number of rows
291 .  n - number of columns
292 .  nz - number of nonzeros per row (same for all rows)
293 -  nnz - array containing the number of nonzeros in the various rows
294          (possibly different for each row) or NULL
295 
296    Output Parameter:
297 .  A - the matrix
298 
299    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(),
300    MatXXXXSetPreallocation() paradigm instead of this routine directly.
301    [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation]
302 
303    Notes:
304    If nnz is given then nz is ignored
305 
306    The AIJ format (also called the Yale sparse matrix format or
307    compressed row storage), is fully compatible with standard Fortran 77
308    storage.  That is, the stored row and column indices can begin at
309    either one (as in Fortran) or zero.  See the users' manual for details.
310 
311    Specify the preallocated storage with either nz or nnz (not both).
312    Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
313    allocation.  For large problems you MUST preallocate memory or you
314    will get TERRIBLE performance, see the users' manual chapter on matrices.
315 
316    Level: intermediate
317 
318 .seealso: MatCreate(), MatCreateAIJ(), MatCreateAIJCUSP(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatCreateAIJ()
319 
320 @*/
321 PetscErrorCode  MatCreateSeqAIJViennaCL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A)
322 {
323   PetscErrorCode ierr;
324 
325   PetscFunctionBegin;
326   ierr = MatCreate(comm,A);CHKERRQ(ierr);
327   ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr);
328   ierr = MatSetType(*A,MATSEQAIJVIENNACL);CHKERRQ(ierr);
329   ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr);
330   PetscFunctionReturn(0);
331 }
332 
333 
334 PetscErrorCode MatDestroy_SeqAIJViennaCL(Mat A)
335 {
336   PetscErrorCode ierr;
337   Mat_SeqAIJViennaCL *viennaclcontainer = (Mat_SeqAIJViennaCL*)A->spptr;
338 
339   PetscFunctionBegin;
340   try {
341     if (viennaclcontainer) {
342       delete viennaclcontainer->tempvec;
343       delete viennaclcontainer->mat;
344       delete viennaclcontainer->compressed_mat;
345       delete viennaclcontainer;
346     }
347     A->valid_GPU_matrix = PETSC_VIENNACL_UNALLOCATED;
348   } catch(std::exception const & ex) {
349     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex.what());
350   }
351 
352   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijviennacl_C",NULL);CHKERRQ(ierr);
353 
354   /* this next line is because MatDestroy tries to PetscFree spptr if it is not zero, and PetscFree only works if the memory was allocated with PetscNew or PetscMalloc, which don't call the constructor */
355   A->spptr = 0;
356   ierr     = MatDestroy_SeqAIJ(A);CHKERRQ(ierr);
357   PetscFunctionReturn(0);
358 }
359 
360 
361 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJViennaCL(Mat B)
362 {
363   PetscErrorCode ierr;
364 
365   PetscFunctionBegin;
366   ierr = MatCreate_SeqAIJ(B);CHKERRQ(ierr);
367   ierr = MatConvert_SeqAIJ_SeqAIJViennaCL(B,MATSEQAIJVIENNACL,MAT_INPLACE_MATRIX,&B);
368   PetscFunctionReturn(0);
369 }
370 
371 static PetscErrorCode MatDuplicate_SeqAIJViennaCL(Mat A,MatDuplicateOption cpvalues,Mat *B)
372 {
373   PetscErrorCode ierr;
374   Mat C;
375 
376   PetscFunctionBegin;
377   ierr = MatDuplicate_SeqAIJ(A,cpvalues,B);CHKERRQ(ierr);
378   C = *B;
379 
380   C->ops->mult        = MatMult_SeqAIJViennaCL;
381   C->ops->multadd     = MatMultAdd_SeqAIJViennaCL;
382   C->ops->assemblyend = MatAssemblyEnd_SeqAIJViennaCL;
383   C->ops->destroy     = MatDestroy_SeqAIJViennaCL;
384   C->ops->getvecs     = MatCreateVecs_SeqAIJViennaCL;
385   C->ops->duplicate   = MatDuplicate_SeqAIJViennaCL;
386 
387   C->spptr        = new Mat_SeqAIJViennaCL();
388   ((Mat_SeqAIJViennaCL*)C->spptr)->tempvec        = NULL;
389   ((Mat_SeqAIJViennaCL*)C->spptr)->mat            = NULL;
390   ((Mat_SeqAIJViennaCL*)C->spptr)->compressed_mat = NULL;
391 
392   ierr = PetscObjectChangeTypeName((PetscObject)C,MATSEQAIJVIENNACL);CHKERRQ(ierr);
393 
394   C->valid_GPU_matrix = PETSC_VIENNACL_UNALLOCATED;
395 
396   /* If the source matrix is already assembled, copy the destination matrix to the GPU */
397   if (C->assembled) {
398     ierr = MatViennaCLCopyToGPU(C);CHKERRQ(ierr);
399   }
400 
401   PetscFunctionReturn(0);
402 }
403 
404 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJViennaCL(Mat A,MatType type,MatReuse reuse,Mat *newmat)
405 {
406   PetscErrorCode ierr;
407   Mat            B;
408   Mat_SeqAIJ     *aij;
409 
410   PetscFunctionBegin;
411 
412   if (reuse == MAT_REUSE_MATRIX) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"MAT_REUSE_MATRIX is not supported. Consider using MAT_INPLACE_MATRIX instead");
413 
414   if (reuse == MAT_INITIAL_MATRIX) {
415     ierr = MatDuplicate(A,MAT_COPY_VALUES,newmat);CHKERRQ(ierr);
416   }
417 
418   B = *newmat;
419 
420   aij             = (Mat_SeqAIJ*)B->data;
421   aij->inode.use  = PETSC_FALSE;
422 
423   B->ops->mult    = MatMult_SeqAIJViennaCL;
424   B->ops->multadd = MatMultAdd_SeqAIJViennaCL;
425   B->spptr        = new Mat_SeqAIJViennaCL();
426 
427   ((Mat_SeqAIJViennaCL*)B->spptr)->tempvec        = NULL;
428   ((Mat_SeqAIJViennaCL*)B->spptr)->mat            = NULL;
429   ((Mat_SeqAIJViennaCL*)B->spptr)->compressed_mat = NULL;
430 
431   B->ops->assemblyend    = MatAssemblyEnd_SeqAIJViennaCL;
432   B->ops->destroy        = MatDestroy_SeqAIJViennaCL;
433   B->ops->getvecs        = MatCreateVecs_SeqAIJViennaCL;
434   B->ops->duplicate      = MatDuplicate_SeqAIJViennaCL;
435 
436   ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJVIENNACL);CHKERRQ(ierr);
437 
438   ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijviennacl_C",MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr);
439 
440   B->valid_GPU_matrix = PETSC_VIENNACL_UNALLOCATED;
441 
442   /* If the source matrix is already assembled, copy the destination matrix to the GPU */
443   if (B->assembled) {
444     ierr = MatViennaCLCopyToGPU(B);CHKERRQ(ierr);
445   }
446 
447   PetscFunctionReturn(0);
448 }
449 
450 
451 /*M
452    MATSEQAIJVIENNACL - MATAIJVIENNACL = "aijviennacl" = "seqaijviennacl" - A matrix type to be used for sparse matrices.
453 
454    A matrix type type whose data resides on GPUs. These matrices are in CSR format by
455    default. All matrix calculations are performed using the ViennaCL library.
456 
457    Options Database Keys:
458 +  -mat_type aijviennacl - sets the matrix type to "seqaijviennacl" during a call to MatSetFromOptions()
459 .  -mat_viennacl_storage_format csr - sets the storage format of matrices for MatMult during a call to MatSetFromOptions().
460 -  -mat_viennacl_mult_storage_format csr - sets the storage format of matrices for MatMult during a call to MatSetFromOptions().
461 
462   Level: beginner
463 
464 .seealso: MatCreateSeqAIJViennaCL(), MATAIJVIENNACL, MatCreateAIJViennaCL()
465 M*/
466 
467 
468 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_ViennaCL(void)
469 {
470   PetscErrorCode ierr;
471 
472   PetscFunctionBegin;
473   ierr = MatSolverPackageRegister(MATSOLVERPETSC, MATSEQAIJVIENNACL,    MAT_FACTOR_LU,MatGetFactor_seqaij_petsc);CHKERRQ(ierr);
474   ierr = MatSolverPackageRegister(MATSOLVERPETSC, MATSEQAIJVIENNACL,    MAT_FACTOR_CHOLESKY,MatGetFactor_seqaij_petsc);CHKERRQ(ierr);
475   ierr = MatSolverPackageRegister(MATSOLVERPETSC, MATSEQAIJVIENNACL,    MAT_FACTOR_ILU,MatGetFactor_seqaij_petsc);CHKERRQ(ierr);
476   ierr = MatSolverPackageRegister(MATSOLVERPETSC, MATSEQAIJVIENNACL,    MAT_FACTOR_ICC,MatGetFactor_seqaij_petsc);CHKERRQ(ierr);
477   PetscFunctionReturn(0);
478 }
479 
480