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