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