1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the LUSOL package of .... 5 6 */ 7 #include "src/mat/impls/aij/seq/aij.h" 8 9 #if defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 10 #define LU1FAC lu1fac_ 11 #define LU6SOL lu6sol_ 12 #define M1PAGE m1page_ 13 #define M5SETX m5setx_ 14 #define M6RDEL m6rdel_ 15 #elif !defined(PETSC_HAVE_FORTRAN_CAPS) 16 #define LU1FAC lu1fac 17 #define LU6SOL lu6sol 18 #define M1PAGE m1page 19 #define M5SETX m5setx 20 #define M6RDEL m6rdel 21 #endif 22 23 EXTERN_C_BEGIN 24 /* 25 Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require 26 */ 27 void PETSC_STDCALL M1PAGE() { 28 ; 29 } 30 void PETSC_STDCALL M5SETX() { 31 ; 32 } 33 34 void PETSC_STDCALL M6RDEL() { 35 ; 36 } 37 38 extern void PETSC_STDCALL LU1FAC (int *m, int *n, int *nnz, int *size, int *luparm, 39 double *parmlu, double *data, int *indc, int *indr, 40 int *rowperm, int *colperm, int *collen, int *rowlen, 41 int *colstart, int *rowstart, int *rploc, int *cploc, 42 int *rpinv, int *cpinv, double *w, int *inform); 43 44 extern void PETSC_STDCALL LU6SOL (int *mode, int *m, int *n, double *rhs, double *x, 45 int *size, int *luparm, double *parmlu, double *data, 46 int *indc, int *indr, int *rowperm, int *colperm, 47 int *collen, int *rowlen, int *colstart, int *rowstart, 48 int *inform); 49 EXTERN_C_END 50 51 EXTERN PetscErrorCode MatDuplicate_LUSOL(Mat,MatDuplicateOption,Mat*); 52 53 typedef struct { 54 double *data; 55 int *indc; 56 int *indr; 57 58 int *ip; 59 int *iq; 60 int *lenc; 61 int *lenr; 62 int *locc; 63 int *locr; 64 int *iploc; 65 int *iqloc; 66 int *ipinv; 67 int *iqinv; 68 double *mnsw; 69 double *mnsv; 70 71 double elbowroom; 72 double luroom; /* Extra space allocated when factor fails */ 73 double parmlu[30]; /* Input/output to LUSOL */ 74 75 int n; /* Number of rows/columns in matrix */ 76 int nz; /* Number of nonzeros */ 77 int nnz; /* Number of nonzeros allocated for factors */ 78 int luparm[30]; /* Input/output to LUSOL */ 79 80 PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*); 81 PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 82 PetscErrorCode (*MatDestroy)(Mat); 83 PetscTruth CleanUpLUSOL; 84 85 } Mat_LUSOL; 86 87 /* LUSOL input/Output Parameters (Description uses C-style indexes 88 * 89 * Input parameters Typical value 90 * 91 * luparm(0) = nout File number for printed messages. 6 92 * luparm(1) = lprint Print level. 0 93 * < 0 suppresses output. 94 * = 0 gives error messages. 95 * = 1 gives debug output from some of the 96 * other routines in LUSOL. 97 * >= 2 gives the pivot row and column and the 98 * no. of rows and columns involved at 99 * each elimination step in lu1fac. 100 * luparm(2) = maxcol lu1fac: maximum number of columns 5 101 * searched allowed in a Markowitz-type 102 * search for the next pivot element. 103 * For some of the factorization, the 104 * number of rows searched is 105 * maxrow = maxcol - 1. 106 * 107 * 108 * Output parameters 109 * 110 * luparm(9) = inform Return code from last call to any LU routine. 111 * luparm(10) = nsing No. of singularities marked in the 112 * output array w(*). 113 * luparm(11) = jsing Column index of last singularity. 114 * luparm(12) = minlen Minimum recommended value for lena. 115 * luparm(13) = maxlen ? 116 * luparm(14) = nupdat No. of updates performed by the lu8 routines. 117 * luparm(15) = nrank No. of nonempty rows of U. 118 * luparm(16) = ndens1 No. of columns remaining when the density of 119 * the matrix being factorized reached dens1. 120 * luparm(17) = ndens2 No. of columns remaining when the density of 121 * the matrix being factorized reached dens2. 122 * luparm(18) = jumin The column index associated with dumin. 123 * luparm(19) = numl0 No. of columns in initial L. 124 * luparm(20) = lenl0 Size of initial L (no. of nonzeros). 125 * luparm(21) = lenu0 Size of initial U. 126 * luparm(22) = lenl Size of current L. 127 * luparm(23) = lenu Size of current U. 128 * luparm(24) = lrow Length of row file. 129 * luparm(25) = ncp No. of compressions of LU data structures. 130 * luparm(26) = mersum lu1fac: sum of Markowitz merit counts. 131 * luparm(27) = nutri lu1fac: triangular rows in U. 132 * luparm(28) = nltri lu1fac: triangular rows in L. 133 * luparm(29) = 134 * 135 * 136 * Input parameters Typical value 137 * 138 * parmlu(0) = elmax1 Max multiplier allowed in L 10.0 139 * during factor. 140 * parmlu(1) = elmax2 Max multiplier allowed in L 10.0 141 * during updates. 142 * parmlu(2) = small Absolute tolerance for eps**0.8 143 * treating reals as zero. IBM double: 3.0d-13 144 * parmlu(3) = utol1 Absolute tol for flagging eps**0.66667 145 * small diagonals of U. IBM double: 3.7d-11 146 * parmlu(4) = utol2 Relative tol for flagging eps**0.66667 147 * small diagonals of U. IBM double: 3.7d-11 148 * parmlu(5) = uspace Factor limiting waste space in U. 3.0 149 * In lu1fac, the row or column lists 150 * are compressed if their length 151 * exceeds uspace times the length of 152 * either file after the last compression. 153 * parmlu(6) = dens1 The density at which the Markowitz 0.3 154 * strategy should search maxcol columns 155 * and no rows. 156 * parmlu(7) = dens2 the density at which the Markowitz 0.6 157 * strategy should search only 1 column 158 * or (preferably) use a dense LU for 159 * all the remaining rows and columns. 160 * 161 * 162 * Output parameters 163 * 164 * parmlu(9) = amax Maximum element in A. 165 * parmlu(10) = elmax Maximum multiplier in current L. 166 * parmlu(11) = umax Maximum element in current U. 167 * parmlu(12) = dumax Maximum diagonal in U. 168 * parmlu(13) = dumin Minimum diagonal in U. 169 * parmlu(14) = 170 * parmlu(15) = 171 * parmlu(16) = 172 * parmlu(17) = 173 * parmlu(18) = 174 * parmlu(19) = resid lu6sol: residual after solve with U or U'. 175 * ... 176 * parmlu(29) = 177 */ 178 179 #define Factorization_Tolerance 1e-1 180 #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0) 181 #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */ 182 183 EXTERN_C_BEGIN 184 #undef __FUNCT__ 185 #define __FUNCT__ "MatConvert_LUSOL_SeqAIJ" 186 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_LUSOL_SeqAIJ(Mat A,const MatType type,MatReuse reuse,Mat *newmat) 187 { 188 PetscErrorCode ierr; 189 Mat B=*newmat; 190 Mat_LUSOL *lusol=(Mat_LUSOL *)A->spptr; 191 192 PetscFunctionBegin; 193 if (reuse == MAT_INITIAL_MATRIX) { 194 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 195 } 196 B->ops->duplicate = lusol->MatDuplicate; 197 B->ops->lufactorsymbolic = lusol->MatLUFactorSymbolic; 198 B->ops->destroy = lusol->MatDestroy; 199 200 ierr = PetscFree(lusol);CHKERRQ(ierr); 201 202 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_lusol_C","",PETSC_NULL);CHKERRQ(ierr); 203 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_lusol_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 204 205 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 206 *newmat = B; 207 PetscFunctionReturn(0); 208 } 209 EXTERN_C_END 210 211 #undef __FUNCT__ 212 #define __FUNCT__ "MatDestroy_LUSOL" 213 PetscErrorCode MatDestroy_LUSOL(Mat A) 214 { 215 PetscErrorCode ierr; 216 Mat_LUSOL *lusol=(Mat_LUSOL *)A->spptr; 217 218 PetscFunctionBegin; 219 if (lusol->CleanUpLUSOL) { 220 ierr = PetscFree(lusol->ip);CHKERRQ(ierr); 221 ierr = PetscFree(lusol->iq);CHKERRQ(ierr); 222 ierr = PetscFree(lusol->lenc);CHKERRQ(ierr); 223 ierr = PetscFree(lusol->lenr);CHKERRQ(ierr); 224 ierr = PetscFree(lusol->locc);CHKERRQ(ierr); 225 ierr = PetscFree(lusol->locr);CHKERRQ(ierr); 226 ierr = PetscFree(lusol->iploc);CHKERRQ(ierr); 227 ierr = PetscFree(lusol->iqloc);CHKERRQ(ierr); 228 ierr = PetscFree(lusol->ipinv);CHKERRQ(ierr); 229 ierr = PetscFree(lusol->iqinv);CHKERRQ(ierr); 230 ierr = PetscFree(lusol->mnsw);CHKERRQ(ierr); 231 ierr = PetscFree(lusol->mnsv);CHKERRQ(ierr); 232 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 233 } 234 235 ierr = MatConvert_LUSOL_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A); 236 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNCT__ 241 #define __FUNCT__ "MatSolve_LUSOL" 242 PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x) 243 { 244 Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr; 245 double *bb,*xx; 246 int mode=5; 247 PetscErrorCode ierr; 248 int i,m,n,nnz,status; 249 250 PetscFunctionBegin; 251 ierr = VecGetArray(x, &xx);CHKERRQ(ierr); 252 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 253 254 m = n = lusol->n; 255 nnz = lusol->nnz; 256 257 for (i = 0; i < m; i++) 258 { 259 lusol->mnsv[i] = bb[i]; 260 } 261 262 LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, 263 lusol->luparm, lusol->parmlu, lusol->data, 264 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 265 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status); 266 267 if (status != 0) 268 { 269 SETERRQ(PETSC_ERR_ARG_SIZ,"solve failed"); 270 } 271 272 ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr); 273 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 274 PetscFunctionReturn(0); 275 } 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "MatLUFactorNumeric_LUSOL" 279 PetscErrorCode MatLUFactorNumeric_LUSOL(Mat A,MatFactorInfo *info,Mat *F) 280 { 281 Mat_SeqAIJ *a; 282 Mat_LUSOL *lusol = (Mat_LUSOL*)(*F)->spptr; 283 PetscErrorCode ierr; 284 int m, n, nz, nnz, status; 285 int i, rs, re; 286 int factorizations; 287 288 PetscFunctionBegin; 289 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr); 290 a = (Mat_SeqAIJ *)A->data; 291 292 if (m != lusol->n) { 293 SETERRQ(PETSC_ERR_ARG_SIZ,"factorization struct inconsistent"); 294 } 295 296 factorizations = 0; 297 do 298 { 299 /*******************************************************************/ 300 /* Check the workspace allocation. */ 301 /*******************************************************************/ 302 303 nz = a->nz; 304 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 305 nnz = PetscMax(nnz, 5*n); 306 307 if (nnz < lusol->luparm[12]){ 308 nnz = (int)(lusol->luroom * lusol->luparm[12]); 309 } else if ((factorizations > 0) && (lusol->luroom < 6)){ 310 lusol->luroom += 0.1; 311 } 312 313 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 314 315 if (nnz > lusol->nnz){ 316 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 317 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);CHKERRQ(ierr); 318 lusol->indr = lusol->indc + nnz; 319 lusol->data = (double *)(lusol->indr + nnz); 320 lusol->nnz = nnz; 321 } 322 323 /*******************************************************************/ 324 /* Fill in the data for the problem. (1-based Fortran style) */ 325 /*******************************************************************/ 326 327 nz = 0; 328 for (i = 0; i < n; i++) 329 { 330 rs = a->i[i]; 331 re = a->i[i+1]; 332 333 while (rs < re) 334 { 335 if (a->a[rs] != 0.0) 336 { 337 lusol->indc[nz] = i + 1; 338 lusol->indr[nz] = a->j[rs] + 1; 339 lusol->data[nz] = a->a[rs]; 340 nz++; 341 } 342 rs++; 343 } 344 } 345 346 /*******************************************************************/ 347 /* Do the factorization. */ 348 /*******************************************************************/ 349 350 LU1FAC(&m, &n, &nz, &nnz, 351 lusol->luparm, lusol->parmlu, lusol->data, 352 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 353 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 354 lusol->iploc, lusol->iqloc, lusol->ipinv, 355 lusol->iqinv, lusol->mnsw, &status); 356 357 switch(status) 358 { 359 case 0: /* factored */ 360 break; 361 362 case 7: /* insufficient memory */ 363 break; 364 365 case 1: 366 case -1: /* singular */ 367 SETERRQ(PETSC_ERR_LIB,"Singular matrix"); 368 369 case 3: 370 case 4: /* error conditions */ 371 SETERRQ(PETSC_ERR_LIB,"matrix error"); 372 373 default: /* unknown condition */ 374 SETERRQ(PETSC_ERR_LIB,"matrix unknown return code"); 375 } 376 377 factorizations++; 378 } while (status == 7); 379 (*F)->assembled = PETSC_TRUE; 380 PetscFunctionReturn(0); 381 } 382 383 #undef __FUNCT__ 384 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL" 385 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) { 386 /************************************************************************/ 387 /* Input */ 388 /* A - matrix to factor */ 389 /* r - row permutation (ignored) */ 390 /* c - column permutation (ignored) */ 391 /* */ 392 /* Output */ 393 /* F - matrix storing the factorization; */ 394 /************************************************************************/ 395 Mat B; 396 Mat_LUSOL *lusol; 397 PetscErrorCode ierr; 398 int i, m, n, nz, nnz; 399 400 PetscFunctionBegin; 401 402 /************************************************************************/ 403 /* Check the arguments. */ 404 /************************************************************************/ 405 406 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 407 nz = ((Mat_SeqAIJ *)A->data)->nz; 408 409 /************************************************************************/ 410 /* Create the factorization. */ 411 /************************************************************************/ 412 413 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 414 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 415 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 416 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 417 418 B->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 419 B->ops->solve = MatSolve_LUSOL; 420 B->factor = FACTOR_LU; 421 lusol = (Mat_LUSOL*)(B->spptr); 422 423 /************************************************************************/ 424 /* Initialize parameters */ 425 /************************************************************************/ 426 427 for (i = 0; i < 30; i++) 428 { 429 lusol->luparm[i] = 0; 430 lusol->parmlu[i] = 0; 431 } 432 433 lusol->luparm[1] = -1; 434 lusol->luparm[2] = 5; 435 lusol->luparm[7] = 1; 436 437 lusol->parmlu[0] = 1 / Factorization_Tolerance; 438 lusol->parmlu[1] = 1 / Factorization_Tolerance; 439 lusol->parmlu[2] = Factorization_Small_Tolerance; 440 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 441 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 442 lusol->parmlu[5] = 3.0; 443 lusol->parmlu[6] = 0.3; 444 lusol->parmlu[7] = 0.6; 445 446 /************************************************************************/ 447 /* Allocate the workspace needed by LUSOL. */ 448 /************************************************************************/ 449 450 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 451 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 452 453 lusol->n = n; 454 lusol->nz = nz; 455 lusol->nnz = nnz; 456 lusol->luroom = 1.75; 457 458 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 459 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 460 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 461 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 462 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 463 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 464 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 465 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 466 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 467 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 468 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 469 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 470 471 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc); 472 lusol->indr = lusol->indc + nnz; 473 lusol->data = (double *)(lusol->indr + nnz); 474 lusol->CleanUpLUSOL = PETSC_TRUE; 475 *F = B; 476 PetscFunctionReturn(0); 477 } 478 479 EXTERN_C_BEGIN 480 #undef __FUNCT__ 481 #define __FUNCT__ "MatConvert_SeqAIJ_LUSOL" 482 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_LUSOL(Mat A,const MatType type,MatReuse reuse,Mat *newmat) 483 { 484 PetscErrorCode ierr; 485 PetscInt m, n; 486 Mat_LUSOL *lusol; 487 Mat B=*newmat; 488 489 PetscFunctionBegin; 490 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 491 if (m != n) { 492 SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 493 } 494 if (reuse == MAT_INITIAL_MATRIX) { 495 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 496 } 497 498 ierr = PetscNewLog(B,Mat_LUSOL,&lusol);CHKERRQ(ierr); 499 lusol->MatDuplicate = A->ops->duplicate; 500 lusol->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 501 lusol->MatDestroy = A->ops->destroy; 502 lusol->CleanUpLUSOL = PETSC_FALSE; 503 504 B->spptr = (void*)lusol; 505 B->ops->duplicate = MatDuplicate_LUSOL; 506 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 507 B->ops->destroy = MatDestroy_LUSOL; 508 509 ierr = PetscInfo(A,"Using LUSOL for LU factorization and solves.\n");CHKERRQ(ierr); 510 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_lusol_C", 511 "MatConvert_SeqAIJ_LUSOL",MatConvert_SeqAIJ_LUSOL);CHKERRQ(ierr); 512 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_lusol_seqaij_C", 513 "MatConvert_LUSOL_SeqAIJ",MatConvert_LUSOL_SeqAIJ);CHKERRQ(ierr); 514 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 515 *newmat = B; 516 PetscFunctionReturn(0); 517 } 518 EXTERN_C_END 519 520 #undef __FUNCT__ 521 #define __FUNCT__ "MatDuplicate_LUSOL" 522 PetscErrorCode MatDuplicate_LUSOL(Mat A, MatDuplicateOption op, Mat *M) { 523 PetscErrorCode ierr; 524 Mat_LUSOL *lu=(Mat_LUSOL *)A->spptr; 525 PetscFunctionBegin; 526 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 527 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_LUSOL));CHKERRQ(ierr); 528 PetscFunctionReturn(0); 529 } 530 531 /*MC 532 MATLUSOL - MATLUSOL = "lusol" - A matrix type providing direct solvers (LU) for sequential matrices 533 via the external package LUSOL. 534 535 If LUSOL is installed (see the manual for 536 instructions on how to declare the existence of external packages), 537 a matrix type can be constructed which invokes LUSOL solvers. 538 After calling MatCreate(...,A), simply call MatSetType(A,MATLUSOL). 539 This matrix type is only supported for double precision real. 540 541 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 542 supported for this matrix type. MatConvert can be called for a fast inplace conversion 543 to and from the MATSEQAIJ matrix type. 544 545 Options Database Keys: 546 . -mat_type lusol - sets the matrix type to "lusol" during a call to MatSetFromOptions() 547 548 Level: beginner 549 550 .seealso: PCLU 551 M*/ 552 553 EXTERN_C_BEGIN 554 #undef __FUNCT__ 555 #define __FUNCT__ "MatCreate_LUSOL" 556 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_LUSOL(Mat A) 557 { 558 PetscErrorCode ierr; 559 560 PetscFunctionBegin; 561 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 562 ierr = MatConvert_SeqAIJ_LUSOL(A,MATLUSOL,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); 563 PetscFunctionReturn(0); 564 } 565 EXTERN_C_END 566