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 /* This routine is only called to convert an unfactored PETSc-LUSOL matrix */ 189 /* to its base PETSc type, so we will ignore 'MatType type'. */ 190 PetscErrorCode ierr; 191 Mat B=*newmat; 192 Mat_LUSOL *lusol=(Mat_LUSOL *)A->spptr; 193 194 PetscFunctionBegin; 195 if (reuse == MAT_INITIAL_MATRIX) { 196 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 197 } 198 B->ops->duplicate = lusol->MatDuplicate; 199 B->ops->lufactorsymbolic = lusol->MatLUFactorSymbolic; 200 B->ops->destroy = lusol->MatDestroy; 201 202 ierr = PetscFree(lusol);CHKERRQ(ierr); 203 204 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_lusol_C","",PETSC_NULL);CHKERRQ(ierr); 205 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_lusol_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 206 207 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 208 *newmat = B; 209 PetscFunctionReturn(0); 210 } 211 EXTERN_C_END 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "MatDestroy_LUSOL" 215 PetscErrorCode MatDestroy_LUSOL(Mat A) 216 { 217 PetscErrorCode ierr; 218 Mat_LUSOL *lusol=(Mat_LUSOL *)A->spptr; 219 220 PetscFunctionBegin; 221 if (lusol->CleanUpLUSOL) { 222 ierr = PetscFree(lusol->ip);CHKERRQ(ierr); 223 ierr = PetscFree(lusol->iq);CHKERRQ(ierr); 224 ierr = PetscFree(lusol->lenc);CHKERRQ(ierr); 225 ierr = PetscFree(lusol->lenr);CHKERRQ(ierr); 226 ierr = PetscFree(lusol->locc);CHKERRQ(ierr); 227 ierr = PetscFree(lusol->locr);CHKERRQ(ierr); 228 ierr = PetscFree(lusol->iploc);CHKERRQ(ierr); 229 ierr = PetscFree(lusol->iqloc);CHKERRQ(ierr); 230 ierr = PetscFree(lusol->ipinv);CHKERRQ(ierr); 231 ierr = PetscFree(lusol->iqinv);CHKERRQ(ierr); 232 ierr = PetscFree(lusol->mnsw);CHKERRQ(ierr); 233 ierr = PetscFree(lusol->mnsv);CHKERRQ(ierr); 234 235 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 236 } 237 238 ierr = MatConvert_LUSOL_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A); 239 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "MatSolve_LUSOL" 245 PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x) 246 { 247 Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr; 248 double *bb,*xx; 249 int mode=5; 250 PetscErrorCode ierr; 251 int i,m,n,nnz,status; 252 253 PetscFunctionBegin; 254 ierr = VecGetArray(x, &xx);CHKERRQ(ierr); 255 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 256 257 m = n = lusol->n; 258 nnz = lusol->nnz; 259 260 for (i = 0; i < m; i++) 261 { 262 lusol->mnsv[i] = bb[i]; 263 } 264 265 LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, 266 lusol->luparm, lusol->parmlu, lusol->data, 267 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 268 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status); 269 270 if (status != 0) 271 { 272 SETERRQ(PETSC_ERR_ARG_SIZ,"solve failed"); 273 } 274 275 ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr); 276 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 277 PetscFunctionReturn(0); 278 } 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "MatLUFactorNumeric_LUSOL" 282 PetscErrorCode MatLUFactorNumeric_LUSOL(Mat A,MatFactorInfo *info,Mat *F) 283 { 284 Mat_SeqAIJ *a; 285 Mat_LUSOL *lusol = (Mat_LUSOL*)(*F)->spptr; 286 PetscErrorCode ierr; 287 int m, n, nz, nnz, status; 288 int i, rs, re; 289 int factorizations; 290 291 PetscFunctionBegin; 292 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr); 293 a = (Mat_SeqAIJ *)A->data; 294 295 if (m != lusol->n) { 296 SETERRQ(PETSC_ERR_ARG_SIZ,"factorization struct inconsistent"); 297 } 298 299 factorizations = 0; 300 do 301 { 302 /*******************************************************************/ 303 /* Check the workspace allocation. */ 304 /*******************************************************************/ 305 306 nz = a->nz; 307 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz)); 308 nnz = PetscMax(nnz, 5*n); 309 310 if (nnz < lusol->luparm[12]){ 311 nnz = (int)(lusol->luroom * lusol->luparm[12]); 312 } else if ((factorizations > 0) && (lusol->luroom < 6)){ 313 lusol->luroom += 0.1; 314 } 315 316 nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23]))); 317 318 if (nnz > lusol->nnz){ 319 ierr = PetscFree(lusol->indc);CHKERRQ(ierr); 320 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);CHKERRQ(ierr); 321 lusol->indr = lusol->indc + nnz; 322 lusol->data = (double *)(lusol->indr + nnz); 323 lusol->nnz = nnz; 324 } 325 326 /*******************************************************************/ 327 /* Fill in the data for the problem. (1-based Fortran style) */ 328 /*******************************************************************/ 329 330 nz = 0; 331 for (i = 0; i < n; i++) 332 { 333 rs = a->i[i]; 334 re = a->i[i+1]; 335 336 while (rs < re) 337 { 338 if (a->a[rs] != 0.0) 339 { 340 lusol->indc[nz] = i + 1; 341 lusol->indr[nz] = a->j[rs] + 1; 342 lusol->data[nz] = a->a[rs]; 343 nz++; 344 } 345 rs++; 346 } 347 } 348 349 /*******************************************************************/ 350 /* Do the factorization. */ 351 /*******************************************************************/ 352 353 LU1FAC(&m, &n, &nz, &nnz, 354 lusol->luparm, lusol->parmlu, lusol->data, 355 lusol->indc, lusol->indr, lusol->ip, lusol->iq, 356 lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, 357 lusol->iploc, lusol->iqloc, lusol->ipinv, 358 lusol->iqinv, lusol->mnsw, &status); 359 360 switch(status) 361 { 362 case 0: /* factored */ 363 break; 364 365 case 7: /* insufficient memory */ 366 break; 367 368 case 1: 369 case -1: /* singular */ 370 SETERRQ(PETSC_ERR_LIB,"Singular matrix"); 371 372 case 3: 373 case 4: /* error conditions */ 374 SETERRQ(PETSC_ERR_LIB,"matrix error"); 375 376 default: /* unknown condition */ 377 SETERRQ(PETSC_ERR_LIB,"matrix unknown return code"); 378 } 379 380 factorizations++; 381 } while (status == 7); 382 (*F)->assembled = PETSC_TRUE; 383 PetscFunctionReturn(0); 384 } 385 386 #undef __FUNCT__ 387 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL" 388 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) { 389 /************************************************************************/ 390 /* Input */ 391 /* A - matrix to factor */ 392 /* r - row permutation (ignored) */ 393 /* c - column permutation (ignored) */ 394 /* */ 395 /* Output */ 396 /* F - matrix storing the factorization; */ 397 /************************************************************************/ 398 Mat B; 399 Mat_LUSOL *lusol; 400 PetscErrorCode ierr; 401 int i, m, n, nz, nnz; 402 403 PetscFunctionBegin; 404 405 /************************************************************************/ 406 /* Check the arguments. */ 407 /************************************************************************/ 408 409 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 410 nz = ((Mat_SeqAIJ *)A->data)->nz; 411 412 /************************************************************************/ 413 /* Create the factorization. */ 414 /************************************************************************/ 415 416 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,m,n,&B);CHKERRQ(ierr); 417 ierr = MatSetType(B,A->type_name);CHKERRQ(ierr); 418 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 419 420 B->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 421 B->ops->solve = MatSolve_LUSOL; 422 B->factor = FACTOR_LU; 423 lusol = (Mat_LUSOL*)(B->spptr); 424 425 /************************************************************************/ 426 /* Initialize parameters */ 427 /************************************************************************/ 428 429 for (i = 0; i < 30; i++) 430 { 431 lusol->luparm[i] = 0; 432 lusol->parmlu[i] = 0; 433 } 434 435 lusol->luparm[1] = -1; 436 lusol->luparm[2] = 5; 437 lusol->luparm[7] = 1; 438 439 lusol->parmlu[0] = 1 / Factorization_Tolerance; 440 lusol->parmlu[1] = 1 / Factorization_Tolerance; 441 lusol->parmlu[2] = Factorization_Small_Tolerance; 442 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 443 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 444 lusol->parmlu[5] = 3.0; 445 lusol->parmlu[6] = 0.3; 446 lusol->parmlu[7] = 0.6; 447 448 /************************************************************************/ 449 /* Allocate the workspace needed by LUSOL. */ 450 /************************************************************************/ 451 452 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 453 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 454 455 lusol->n = n; 456 lusol->nz = nz; 457 lusol->nnz = nnz; 458 lusol->luroom = 1.75; 459 460 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 461 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 462 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 463 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 464 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 465 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 466 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 467 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 468 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 469 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 470 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 471 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 472 473 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc); 474 lusol->indr = lusol->indc + nnz; 475 lusol->data = (double *)(lusol->indr + nnz); 476 lusol->CleanUpLUSOL = PETSC_TRUE; 477 *F = B; 478 PetscFunctionReturn(0); 479 } 480 481 EXTERN_C_BEGIN 482 #undef __FUNCT__ 483 #define __FUNCT__ "MatConvert_SeqAIJ_LUSOL" 484 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_LUSOL(Mat A,const MatType type,MatReuse reuse,Mat *newmat) 485 { 486 PetscErrorCode ierr; 487 PetscInt m, n; 488 Mat_LUSOL *lusol; 489 Mat B=*newmat; 490 491 PetscFunctionBegin; 492 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 493 if (m != n) { 494 SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 495 } 496 if (reuse == MAT_INITIAL_MATRIX) { 497 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 498 } 499 500 ierr = PetscNew(Mat_LUSOL,&lusol);CHKERRQ(ierr); 501 lusol->MatDuplicate = A->ops->duplicate; 502 lusol->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 503 lusol->MatDestroy = A->ops->destroy; 504 lusol->CleanUpLUSOL = PETSC_FALSE; 505 506 B->spptr = (void*)lusol; 507 B->ops->duplicate = MatDuplicate_LUSOL; 508 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 509 B->ops->destroy = MatDestroy_LUSOL; 510 511 ierr = PetscLogInfo((0,"MatConvert_SeqAIJ_LUSOL:Using LUSOL for LU factorization and solves.\n"));CHKERRQ(ierr); 512 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_lusol_C", 513 "MatConvert_SeqAIJ_LUSOL",MatConvert_SeqAIJ_LUSOL);CHKERRQ(ierr); 514 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_lusol_seqaij_C", 515 "MatConvert_LUSOL_SeqAIJ",MatConvert_LUSOL_SeqAIJ);CHKERRQ(ierr); 516 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 517 *newmat = B; 518 PetscFunctionReturn(0); 519 } 520 EXTERN_C_END 521 522 #undef __FUNCT__ 523 #define __FUNCT__ "MatDuplicate_LUSOL" 524 PetscErrorCode MatDuplicate_LUSOL(Mat A, MatDuplicateOption op, Mat *M) { 525 PetscErrorCode ierr; 526 Mat_LUSOL *lu=(Mat_LUSOL *)A->spptr; 527 PetscFunctionBegin; 528 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 529 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_LUSOL));CHKERRQ(ierr); 530 PetscFunctionReturn(0); 531 } 532 533 /*MC 534 MATLUSOL - MATLUSOL = "lusol" - A matrix type providing direct solvers (LU) for sequential matrices 535 via the external package LUSOL. 536 537 If LUSOL is installed (see the manual for 538 instructions on how to declare the existence of external packages), 539 a matrix type can be constructed which invokes LUSOL solvers. 540 After calling MatCreate(...,A), simply call MatSetType(A,MATLUSOL). 541 This matrix type is only supported for double precision real. 542 543 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 544 supported for this matrix type. MatConvert can be called for a fast inplace conversion 545 to and from the MATSEQAIJ matrix type. 546 547 Options Database Keys: 548 . -mat_type lusol - sets the matrix type to "lusol" during a call to MatSetFromOptions() 549 550 Level: beginner 551 552 .seealso: PCLU 553 M*/ 554 555 EXTERN_C_BEGIN 556 #undef __FUNCT__ 557 #define __FUNCT__ "MatCreate_LUSOL" 558 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_LUSOL(Mat A) 559 { 560 PetscErrorCode ierr; 561 562 PetscFunctionBegin; 563 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and LUSOL types */ 564 ierr = PetscObjectChangeTypeName((PetscObject)A,MATLUSOL);CHKERRQ(ierr); 565 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 566 ierr = MatConvert_SeqAIJ_LUSOL(A,MATLUSOL,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); 567 PetscFunctionReturn(0); 568 } 569 EXTERN_C_END 570