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