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,&B);CHKERRQ(ierr); 417 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 418 ierr = MatSetType(B,A->type_name);CHKERRQ(ierr); 419 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 420 421 B->ops->lufactornumeric = MatLUFactorNumeric_LUSOL; 422 B->ops->solve = MatSolve_LUSOL; 423 B->factor = FACTOR_LU; 424 lusol = (Mat_LUSOL*)(B->spptr); 425 426 /************************************************************************/ 427 /* Initialize parameters */ 428 /************************************************************************/ 429 430 for (i = 0; i < 30; i++) 431 { 432 lusol->luparm[i] = 0; 433 lusol->parmlu[i] = 0; 434 } 435 436 lusol->luparm[1] = -1; 437 lusol->luparm[2] = 5; 438 lusol->luparm[7] = 1; 439 440 lusol->parmlu[0] = 1 / Factorization_Tolerance; 441 lusol->parmlu[1] = 1 / Factorization_Tolerance; 442 lusol->parmlu[2] = Factorization_Small_Tolerance; 443 lusol->parmlu[3] = Factorization_Pivot_Tolerance; 444 lusol->parmlu[4] = Factorization_Pivot_Tolerance; 445 lusol->parmlu[5] = 3.0; 446 lusol->parmlu[6] = 0.3; 447 lusol->parmlu[7] = 0.6; 448 449 /************************************************************************/ 450 /* Allocate the workspace needed by LUSOL. */ 451 /************************************************************************/ 452 453 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill); 454 nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n); 455 456 lusol->n = n; 457 lusol->nz = nz; 458 lusol->nnz = nnz; 459 lusol->luroom = 1.75; 460 461 ierr = PetscMalloc(sizeof(int)*n,&lusol->ip); 462 ierr = PetscMalloc(sizeof(int)*n,&lusol->iq); 463 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenc); 464 ierr = PetscMalloc(sizeof(int)*n,&lusol->lenr); 465 ierr = PetscMalloc(sizeof(int)*n,&lusol->locc); 466 ierr = PetscMalloc(sizeof(int)*n,&lusol->locr); 467 ierr = PetscMalloc(sizeof(int)*n,&lusol->iploc); 468 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqloc); 469 ierr = PetscMalloc(sizeof(int)*n,&lusol->ipinv); 470 ierr = PetscMalloc(sizeof(int)*n,&lusol->iqinv); 471 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsw); 472 ierr = PetscMalloc(sizeof(double)*n,&lusol->mnsv); 473 474 ierr = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc); 475 lusol->indr = lusol->indc + nnz; 476 lusol->data = (double *)(lusol->indr + nnz); 477 lusol->CleanUpLUSOL = PETSC_TRUE; 478 *F = B; 479 PetscFunctionReturn(0); 480 } 481 482 EXTERN_C_BEGIN 483 #undef __FUNCT__ 484 #define __FUNCT__ "MatConvert_SeqAIJ_LUSOL" 485 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_LUSOL(Mat A,const MatType type,MatReuse reuse,Mat *newmat) 486 { 487 PetscErrorCode ierr; 488 PetscInt m, n; 489 Mat_LUSOL *lusol; 490 Mat B=*newmat; 491 492 PetscFunctionBegin; 493 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 494 if (m != n) { 495 SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 496 } 497 if (reuse == MAT_INITIAL_MATRIX) { 498 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 499 } 500 501 ierr = PetscNew(Mat_LUSOL,&lusol);CHKERRQ(ierr); 502 lusol->MatDuplicate = A->ops->duplicate; 503 lusol->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 504 lusol->MatDestroy = A->ops->destroy; 505 lusol->CleanUpLUSOL = PETSC_FALSE; 506 507 B->spptr = (void*)lusol; 508 B->ops->duplicate = MatDuplicate_LUSOL; 509 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL; 510 B->ops->destroy = MatDestroy_LUSOL; 511 512 ierr = PetscLogInfo((0,"MatConvert_SeqAIJ_LUSOL:Using LUSOL for LU factorization and solves.\n"));CHKERRQ(ierr); 513 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_lusol_C", 514 "MatConvert_SeqAIJ_LUSOL",MatConvert_SeqAIJ_LUSOL);CHKERRQ(ierr); 515 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_lusol_seqaij_C", 516 "MatConvert_LUSOL_SeqAIJ",MatConvert_LUSOL_SeqAIJ);CHKERRQ(ierr); 517 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 518 *newmat = B; 519 PetscFunctionReturn(0); 520 } 521 EXTERN_C_END 522 523 #undef __FUNCT__ 524 #define __FUNCT__ "MatDuplicate_LUSOL" 525 PetscErrorCode MatDuplicate_LUSOL(Mat A, MatDuplicateOption op, Mat *M) { 526 PetscErrorCode ierr; 527 Mat_LUSOL *lu=(Mat_LUSOL *)A->spptr; 528 PetscFunctionBegin; 529 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 530 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_LUSOL));CHKERRQ(ierr); 531 PetscFunctionReturn(0); 532 } 533 534 /*MC 535 MATLUSOL - MATLUSOL = "lusol" - A matrix type providing direct solvers (LU) for sequential matrices 536 via the external package LUSOL. 537 538 If LUSOL is installed (see the manual for 539 instructions on how to declare the existence of external packages), 540 a matrix type can be constructed which invokes LUSOL solvers. 541 After calling MatCreate(...,A), simply call MatSetType(A,MATLUSOL). 542 This matrix type is only supported for double precision real. 543 544 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 545 supported for this matrix type. MatConvert can be called for a fast inplace conversion 546 to and from the MATSEQAIJ matrix type. 547 548 Options Database Keys: 549 . -mat_type lusol - sets the matrix type to "lusol" during a call to MatSetFromOptions() 550 551 Level: beginner 552 553 .seealso: PCLU 554 M*/ 555 556 EXTERN_C_BEGIN 557 #undef __FUNCT__ 558 #define __FUNCT__ "MatCreate_LUSOL" 559 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_LUSOL(Mat A) 560 { 561 PetscErrorCode ierr; 562 563 PetscFunctionBegin; 564 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and LUSOL types */ 565 ierr = PetscObjectChangeTypeName((PetscObject)A,MATLUSOL);CHKERRQ(ierr); 566 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 567 ierr = MatConvert_SeqAIJ_LUSOL(A,MATLUSOL,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); 568 PetscFunctionReturn(0); 569 } 570 EXTERN_C_END 571