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