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 *ygpu = viennacl::linalg::prod(*viennaclstruct->mat,*xgpu); 184 ViennaCLWaitForGPU(); 185 } catch (std::exception const & ex) { 186 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex.what()); 187 } 188 ierr = VecViennaCLRestoreArrayRead(xx,&xgpu);CHKERRQ(ierr); 189 ierr = VecViennaCLRestoreArrayWrite(yy,&ygpu);CHKERRQ(ierr); 190 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 191 } else { 192 ierr = VecSet(yy,0);CHKERRQ(ierr); 193 } 194 PetscFunctionReturn(0); 195 } 196 197 PetscErrorCode MatMultAdd_SeqAIJViennaCL(Mat A,Vec xx,Vec yy,Vec zz) 198 { 199 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 200 PetscErrorCode ierr; 201 Mat_SeqAIJViennaCL *viennaclstruct = (Mat_SeqAIJViennaCL*)A->spptr; 202 const ViennaCLVector *xgpu=NULL,*ygpu=NULL; 203 ViennaCLVector *zgpu=NULL; 204 205 PetscFunctionBegin; 206 if (A->rmap->n > 0 && A->cmap->n > 0 && a->nz) { 207 try { 208 ierr = VecViennaCLGetArrayRead(xx,&xgpu);CHKERRQ(ierr); 209 ierr = VecViennaCLGetArrayRead(yy,&ygpu);CHKERRQ(ierr); 210 ierr = VecViennaCLGetArrayWrite(zz,&zgpu);CHKERRQ(ierr); 211 212 if (a->compressedrow.use) { 213 ViennaCLVector temp = viennacl::linalg::prod(*viennaclstruct->compressed_mat, *xgpu); 214 *zgpu = *ygpu + temp; 215 ViennaCLWaitForGPU(); 216 } else { 217 if (zz == xx || zz == yy) { //temporary required 218 ViennaCLVector temp = viennacl::linalg::prod(*viennaclstruct->mat, *xgpu); 219 *zgpu = *ygpu; 220 *zgpu += temp; 221 ViennaCLWaitForGPU(); 222 } else { 223 *viennaclstruct->tempvec = viennacl::linalg::prod(*viennaclstruct->mat, *xgpu); 224 *zgpu = *ygpu + *viennaclstruct->tempvec; 225 ViennaCLWaitForGPU(); 226 } 227 } 228 229 ierr = VecViennaCLRestoreArrayRead(xx,&xgpu);CHKERRQ(ierr); 230 ierr = VecViennaCLRestoreArrayRead(yy,&ygpu);CHKERRQ(ierr); 231 ierr = VecViennaCLRestoreArrayWrite(zz,&zgpu);CHKERRQ(ierr); 232 233 } catch(std::exception const & ex) { 234 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex.what()); 235 } 236 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 237 } else { 238 ierr = VecCopy(yy,zz);CHKERRQ(ierr); 239 } 240 PetscFunctionReturn(0); 241 } 242 243 PetscErrorCode MatAssemblyEnd_SeqAIJViennaCL(Mat A,MatAssemblyType mode) 244 { 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 ierr = MatAssemblyEnd_SeqAIJ(A,mode);CHKERRQ(ierr); 249 ierr = MatViennaCLCopyToGPU(A);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 /* --------------------------------------------------------------------------------*/ 254 /*@ 255 MatCreateSeqAIJViennaCL - Creates a sparse matrix in AIJ (compressed row) format 256 (the default parallel PETSc format). This matrix will ultimately be pushed down 257 to GPUs and use the ViennaCL library for calculations. For good matrix 258 assembly performance the user should preallocate the matrix storage by setting 259 the parameter nz (or the array nnz). By setting these parameters accurately, 260 performance during matrix assembly can be increased substantially. 261 262 263 Collective on MPI_Comm 264 265 Input Parameters: 266 + comm - MPI communicator, set to PETSC_COMM_SELF 267 . m - number of rows 268 . n - number of columns 269 . nz - number of nonzeros per row (same for all rows) 270 - nnz - array containing the number of nonzeros in the various rows 271 (possibly different for each row) or NULL 272 273 Output Parameter: 274 . A - the matrix 275 276 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 277 MatXXXXSetPreallocation() paradigm instead of this routine directly. 278 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 279 280 Notes: 281 If nnz is given then nz is ignored 282 283 The AIJ format (also called the Yale sparse matrix format or 284 compressed row storage), is fully compatible with standard Fortran 77 285 storage. That is, the stored row and column indices can begin at 286 either one (as in Fortran) or zero. See the users' manual for details. 287 288 Specify the preallocated storage with either nz or nnz (not both). 289 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 290 allocation. For large problems you MUST preallocate memory or you 291 will get TERRIBLE performance, see the users' manual chapter on matrices. 292 293 Level: intermediate 294 295 .seealso: MatCreate(), MatCreateAIJ(), MatCreateAIJCUSPARSE(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatCreateAIJ() 296 297 @*/ 298 PetscErrorCode MatCreateSeqAIJViennaCL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 299 { 300 PetscErrorCode ierr; 301 302 PetscFunctionBegin; 303 ierr = MatCreate(comm,A);CHKERRQ(ierr); 304 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 305 ierr = MatSetType(*A,MATSEQAIJVIENNACL);CHKERRQ(ierr); 306 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 307 PetscFunctionReturn(0); 308 } 309 310 311 PetscErrorCode MatDestroy_SeqAIJViennaCL(Mat A) 312 { 313 PetscErrorCode ierr; 314 Mat_SeqAIJViennaCL *viennaclcontainer = (Mat_SeqAIJViennaCL*)A->spptr; 315 316 PetscFunctionBegin; 317 try { 318 if (viennaclcontainer) { 319 delete viennaclcontainer->tempvec; 320 delete viennaclcontainer->mat; 321 delete viennaclcontainer->compressed_mat; 322 delete viennaclcontainer; 323 } 324 A->valid_GPU_matrix = PETSC_OFFLOAD_UNALLOCATED; 325 } catch(std::exception const & ex) { 326 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex.what()); 327 } 328 329 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijviennacl_C",NULL);CHKERRQ(ierr); 330 331 /* 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 */ 332 A->spptr = 0; 333 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 334 PetscFunctionReturn(0); 335 } 336 337 338 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJViennaCL(Mat B) 339 { 340 PetscErrorCode ierr; 341 342 PetscFunctionBegin; 343 ierr = MatCreate_SeqAIJ(B);CHKERRQ(ierr); 344 ierr = MatConvert_SeqAIJ_SeqAIJViennaCL(B,MATSEQAIJVIENNACL,MAT_INPLACE_MATRIX,&B); 345 PetscFunctionReturn(0); 346 } 347 348 static PetscErrorCode MatDuplicate_SeqAIJViennaCL(Mat A,MatDuplicateOption cpvalues,Mat *B) 349 { 350 PetscErrorCode ierr; 351 Mat C; 352 353 PetscFunctionBegin; 354 ierr = MatDuplicate_SeqAIJ(A,cpvalues,B);CHKERRQ(ierr); 355 C = *B; 356 357 C->ops->mult = MatMult_SeqAIJViennaCL; 358 C->ops->multadd = MatMultAdd_SeqAIJViennaCL; 359 C->ops->assemblyend = MatAssemblyEnd_SeqAIJViennaCL; 360 C->ops->destroy = MatDestroy_SeqAIJViennaCL; 361 C->ops->duplicate = MatDuplicate_SeqAIJViennaCL; 362 363 C->spptr = new Mat_SeqAIJViennaCL(); 364 ((Mat_SeqAIJViennaCL*)C->spptr)->tempvec = NULL; 365 ((Mat_SeqAIJViennaCL*)C->spptr)->mat = NULL; 366 ((Mat_SeqAIJViennaCL*)C->spptr)->compressed_mat = NULL; 367 368 ierr = PetscObjectChangeTypeName((PetscObject)C,MATSEQAIJVIENNACL);CHKERRQ(ierr); 369 370 C->valid_GPU_matrix = PETSC_OFFLOAD_UNALLOCATED; 371 372 /* If the source matrix is already assembled, copy the destination matrix to the GPU */ 373 if (C->assembled) { 374 ierr = MatViennaCLCopyToGPU(C);CHKERRQ(ierr); 375 } 376 377 PetscFunctionReturn(0); 378 } 379 380 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJViennaCL(Mat A,MatType type,MatReuse reuse,Mat *newmat) 381 { 382 PetscErrorCode ierr; 383 Mat B; 384 Mat_SeqAIJ *aij; 385 386 PetscFunctionBegin; 387 388 if (reuse == MAT_REUSE_MATRIX) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"MAT_REUSE_MATRIX is not supported. Consider using MAT_INPLACE_MATRIX instead"); 389 390 if (reuse == MAT_INITIAL_MATRIX) { 391 ierr = MatDuplicate(A,MAT_COPY_VALUES,newmat);CHKERRQ(ierr); 392 } 393 394 B = *newmat; 395 396 aij = (Mat_SeqAIJ*)B->data; 397 aij->inode.use = PETSC_FALSE; 398 399 B->ops->mult = MatMult_SeqAIJViennaCL; 400 B->ops->multadd = MatMultAdd_SeqAIJViennaCL; 401 B->spptr = new Mat_SeqAIJViennaCL(); 402 403 ((Mat_SeqAIJViennaCL*)B->spptr)->tempvec = NULL; 404 ((Mat_SeqAIJViennaCL*)B->spptr)->mat = NULL; 405 ((Mat_SeqAIJViennaCL*)B->spptr)->compressed_mat = NULL; 406 407 B->ops->assemblyend = MatAssemblyEnd_SeqAIJViennaCL; 408 B->ops->destroy = MatDestroy_SeqAIJViennaCL; 409 B->ops->duplicate = MatDuplicate_SeqAIJViennaCL; 410 411 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJVIENNACL);CHKERRQ(ierr); 412 ierr = PetscFree(B->defaultvectype);CHKERRQ(ierr); 413 ierr = PetscStrallocpy(VECVIENNACL,&B->defaultvectype);CHKERRQ(ierr); 414 415 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijviennacl_C",MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 416 417 B->valid_GPU_matrix = PETSC_OFFLOAD_UNALLOCATED; 418 419 /* If the source matrix is already assembled, copy the destination matrix to the GPU */ 420 if (B->assembled) { 421 ierr = MatViennaCLCopyToGPU(B);CHKERRQ(ierr); 422 } 423 424 PetscFunctionReturn(0); 425 } 426 427 428 /*MC 429 MATSEQAIJVIENNACL - MATAIJVIENNACL = "aijviennacl" = "seqaijviennacl" - A matrix type to be used for sparse matrices. 430 431 A matrix type type whose data resides on GPUs. These matrices are in CSR format by 432 default. All matrix calculations are performed using the ViennaCL library. 433 434 Options Database Keys: 435 + -mat_type aijviennacl - sets the matrix type to "seqaijviennacl" during a call to MatSetFromOptions() 436 . -mat_viennacl_storage_format csr - sets the storage format of matrices for MatMult during a call to MatSetFromOptions(). 437 - -mat_viennacl_mult_storage_format csr - sets the storage format of matrices for MatMult during a call to MatSetFromOptions(). 438 439 Level: beginner 440 441 .seealso: MatCreateSeqAIJViennaCL(), MATAIJVIENNACL, MatCreateAIJViennaCL() 442 M*/ 443 444 445 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_ViennaCL(void) 446 { 447 PetscErrorCode ierr; 448 449 PetscFunctionBegin; 450 ierr = MatSolverTypeRegister(MATSOLVERPETSC, MATSEQAIJVIENNACL, MAT_FACTOR_LU,MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 451 ierr = MatSolverTypeRegister(MATSOLVERPETSC, MATSEQAIJVIENNACL, MAT_FACTOR_CHOLESKY,MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 452 ierr = MatSolverTypeRegister(MATSOLVERPETSC, MATSEQAIJVIENNACL, MAT_FACTOR_ILU,MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 453 ierr = MatSolverTypeRegister(MATSOLVERPETSC, MATSEQAIJVIENNACL, MAT_FACTOR_ICC,MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 454 PetscFunctionReturn(0); 455 } 456