xref: /petsc/src/mat/impls/aij/seq/lusol/lusol.c (revision ee45ca4afdde1af4d1deda7cd4dc1a4a63a3ea97)
1 /*
2         Provides an interface to the LUSOL package of ....
3 
4 */
5 #include "src/mat/impls/aij/seq/aij.h"
6 
7 #if defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
8 #define LU1FAC   lu1fac_
9 #define LU6SOL   lu6sol_
10 #define M1PAGE   m1page_
11 #define M5SETX   m5setx_
12 #define M6RDEL   m6rdel_
13 #elif !defined(PETSC_HAVE_FORTRAN_CAPS)
14 #define LU1FAC   lu1fac
15 #define LU6SOL   lu6sol
16 #define M1PAGE   m1page
17 #define M5SETX   m5setx
18 #define M6RDEL   m6rdel
19 #endif
20 
21 EXTERN_C_BEGIN
22 /*
23     Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require
24 */
25 void PETSC_STDCALL M1PAGE() {
26   ;
27 }
28 void PETSC_STDCALL M5SETX() {
29   ;
30 }
31 
32 void PETSC_STDCALL M6RDEL() {
33   ;
34 }
35 
36 extern void PETSC_STDCALL LU1FAC (int *m, int *n, int *nnz, int *size, int *luparm,
37                         double *parmlu, double *data, int *indc, int *indr,
38                         int *rowperm, int *colperm, int *collen, int *rowlen,
39                         int *colstart, int *rowstart, int *rploc, int *cploc,
40                         int *rpinv, int *cpinv, double *w, int *inform);
41 
42 extern void PETSC_STDCALL LU6SOL (int *mode, int *m, int *n, double *rhs, double *x,
43                         int *size, int *luparm, double *parmlu, double *data,
44                         int *indc, int *indr, int *rowperm, int *colperm,
45                         int *collen, int *rowlen, int *colstart, int *rowstart,
46                         int *inform);
47 EXTERN_C_END
48 
49 EXTERN PetscErrorCode MatDuplicate_LUSOL(Mat,MatDuplicateOption,Mat*);
50 
51 typedef struct  {
52   double *data;
53   int *indc;
54   int *indr;
55 
56   int *ip;
57   int *iq;
58   int *lenc;
59   int *lenr;
60   int *locc;
61   int *locr;
62   int *iploc;
63   int *iqloc;
64   int *ipinv;
65   int *iqinv;
66   double *mnsw;
67   double *mnsv;
68 
69   double elbowroom;
70   double luroom;		/* Extra space allocated when factor fails   */
71   double parmlu[30];		/* Input/output to LUSOL                     */
72 
73   int n;			/* Number of rows/columns in matrix          */
74   int nz;			/* Number of nonzeros                        */
75   int nnz;			/* Number of nonzeros allocated for factors  */
76   int luparm[30];		/* Input/output to LUSOL                     */
77 
78   PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
79   PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
80   PetscErrorCode (*MatDestroy)(Mat);
81   PetscTruth CleanUpLUSOL;
82 
83 } Mat_LUSOL;
84 
85 /*  LUSOL input/Output Parameters (Description uses C-style indexes
86  *
87  *  Input parameters                                        Typical value
88  *
89  *  luparm(0) = nout     File number for printed messages.         6
90  *  luparm(1) = lprint   Print level.                              0
91  *                    < 0 suppresses output.
92  *                    = 0 gives error messages.
93  *                    = 1 gives debug output from some of the
94  *                        other routines in LUSOL.
95  *                   >= 2 gives the pivot row and column and the
96  *                        no. of rows and columns involved at
97  *                        each elimination step in lu1fac.
98  *  luparm(2) = maxcol   lu1fac: maximum number of columns         5
99  *                        searched allowed in a Markowitz-type
100  *                        search for the next pivot element.
101  *                        For some of the factorization, the
102  *                        number of rows searched is
103  *                        maxrow = maxcol - 1.
104  *
105  *
106  *  Output parameters
107  *
108  *  luparm(9) = inform   Return code from last call to any LU routine.
109  *  luparm(10) = nsing    No. of singularities marked in the
110  *                        output array w(*).
111  *  luparm(11) = jsing    Column index of last singularity.
112  *  luparm(12) = minlen   Minimum recommended value for  lena.
113  *  luparm(13) = maxlen   ?
114  *  luparm(14) = nupdat   No. of updates performed by the lu8 routines.
115  *  luparm(15) = nrank    No. of nonempty rows of U.
116  *  luparm(16) = ndens1   No. of columns remaining when the density of
117  *                        the matrix being factorized reached dens1.
118  *  luparm(17) = ndens2   No. of columns remaining when the density of
119  *                        the matrix being factorized reached dens2.
120  *  luparm(18) = jumin    The column index associated with dumin.
121  *  luparm(19) = numl0    No. of columns in initial  L.
122  *  luparm(20) = lenl0    Size of initial  L  (no. of nonzeros).
123  *  luparm(21) = lenu0    Size of initial  U.
124  *  luparm(22) = lenl     Size of current  L.
125  *  luparm(23) = lenu     Size of current  U.
126  *  luparm(24) = lrow     Length of row file.
127  *  luparm(25) = ncp      No. of compressions of LU data structures.
128  *  luparm(26) = mersum   lu1fac: sum of Markowitz merit counts.
129  *  luparm(27) = nutri    lu1fac: triangular rows in U.
130  *  luparm(28) = nltri    lu1fac: triangular rows in L.
131  *  luparm(29) =
132  *
133  *
134  *  Input parameters                                        Typical value
135  *
136  *  parmlu(0) = elmax1   Max multiplier allowed in  L           10.0
137  *                        during factor.
138  *  parmlu(1) = elmax2   Max multiplier allowed in  L           10.0
139  *                        during updates.
140  *  parmlu(2) = small    Absolute tolerance for             eps**0.8
141  *                        treating reals as zero.     IBM double: 3.0d-13
142  *  parmlu(3) = utol1    Absolute tol for flagging          eps**0.66667
143  *                        small diagonals of U.       IBM double: 3.7d-11
144  *  parmlu(4) = utol2    Relative tol for flagging          eps**0.66667
145  *                        small diagonals of U.       IBM double: 3.7d-11
146  *  parmlu(5) = uspace   Factor limiting waste space in  U.      3.0
147  *                        In lu1fac, the row or column lists
148  *                        are compressed if their length
149  *                        exceeds uspace times the length of
150  *                        either file after the last compression.
151  *  parmlu(6) = dens1    The density at which the Markowitz      0.3
152  *                        strategy should search maxcol columns
153  *                        and no rows.
154  *  parmlu(7) = dens2    the density at which the Markowitz      0.6
155  *                        strategy should search only 1 column
156  *                        or (preferably) use a dense LU for
157  *                        all the remaining rows and columns.
158  *
159  *
160  *  Output parameters
161  *
162  *  parmlu(9) = amax     Maximum element in  A.
163  *  parmlu(10) = elmax    Maximum multiplier in current  L.
164  *  parmlu(11) = umax     Maximum element in current  U.
165  *  parmlu(12) = dumax    Maximum diagonal in  U.
166  *  parmlu(13) = dumin    Minimum diagonal in  U.
167  *  parmlu(14) =
168  *  parmlu(15) =
169  *  parmlu(16) =
170  *  parmlu(17) =
171  *  parmlu(18) =
172  *  parmlu(19) = resid    lu6sol: residual after solve with U or U'.
173  *  ...
174  *  parmlu(29) =
175  */
176 
177 #define Factorization_Tolerance       1e-1
178 #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0)
179 #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */
180 
181 EXTERN_C_BEGIN
182 #undef __FUNCT__
183 #define __FUNCT__ "MatConvert_LUSOL_SeqAIJ"
184 PetscErrorCode MatConvert_LUSOL_SeqAIJ(Mat A,const MatType type,Mat *newmat)
185 {
186   /* This routine is only called to convert an unfactored PETSc-LUSOL matrix */
187   /* to its base PETSc type, so we will ignore 'MatType type'. */
188   PetscErrorCode ierr;
189   Mat            B=*newmat;
190   Mat_LUSOL      *lusol=(Mat_LUSOL *)A->spptr;
191 
192   PetscFunctionBegin;
193   if (B != A) {
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 
233     ierr = PetscFree(lusol->indc);CHKERRQ(ierr);
234   }
235 
236   ierr = MatConvert_LUSOL_SeqAIJ(A,MATSEQAIJ,&A);
237   ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
238   PetscFunctionReturn(0);
239 }
240 
241 #undef __FUNCT__
242 #define __FUNCT__  "MatSolve_LUSOL"
243 PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x)
244 {
245   Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr;
246   double    *bb,*xx;
247   int       mode=5;
248   PetscErrorCode ierr;
249   int       i,m,n,nnz,status;
250 
251   PetscFunctionBegin;
252   ierr = VecGetArray(x, &xx);CHKERRQ(ierr);
253   ierr = VecGetArray(b, &bb);CHKERRQ(ierr);
254 
255   m = n = lusol->n;
256   nnz = lusol->nnz;
257 
258   for (i = 0; i < m; i++)
259     {
260       lusol->mnsv[i] = bb[i];
261     }
262 
263   LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz,
264          lusol->luparm, lusol->parmlu, lusol->data,
265          lusol->indc, lusol->indr, lusol->ip, lusol->iq,
266          lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status);
267 
268   if (status != 0)
269     {
270       SETERRQ(PETSC_ERR_ARG_SIZ,"solve failed");
271     }
272 
273   ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr);
274   ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr);
275   PetscFunctionReturn(0);
276 }
277 
278 #undef __FUNCT__
279 #define __FUNCT__ "MatLUFactorNumeric_LUSOL"
280 PetscErrorCode MatLUFactorNumeric_LUSOL(Mat A,MatFactorInfo *info,Mat *F)
281 {
282   Mat_SeqAIJ     *a;
283   Mat_LUSOL      *lusol = (Mat_LUSOL*)(*F)->spptr;
284   PetscErrorCode ierr;
285   int            m, n, nz, nnz, status;
286   int            i, rs, re;
287   int            factorizations;
288 
289   PetscFunctionBegin;
290   ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr);
291   a = (Mat_SeqAIJ *)A->data;
292 
293   if (m != lusol->n) {
294     SETERRQ(PETSC_ERR_ARG_SIZ,"factorization struct inconsistent");
295   }
296 
297   factorizations = 0;
298   do
299     {
300       /*******************************************************************/
301       /* Check the workspace allocation.                                 */
302       /*******************************************************************/
303 
304       nz = a->nz;
305       nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz));
306       nnz = PetscMax(nnz, 5*n);
307 
308       if (nnz < lusol->luparm[12]){
309         nnz = (int)(lusol->luroom * lusol->luparm[12]);
310       } else if ((factorizations > 0) && (lusol->luroom < 6)){
311         lusol->luroom += 0.1;
312       }
313 
314       nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23])));
315 
316       if (nnz > lusol->nnz){
317         ierr = PetscFree(lusol->indc);CHKERRQ(ierr);
318         ierr        = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);CHKERRQ(ierr);
319         lusol->indr = lusol->indc + nnz;
320         lusol->data = (double *)(lusol->indr + nnz);
321         lusol->nnz  = nnz;
322       }
323 
324       /*******************************************************************/
325       /* Fill in the data for the problem.      (1-based Fortran style)  */
326       /*******************************************************************/
327 
328       nz = 0;
329       for (i = 0; i < n; i++)
330         {
331           rs = a->i[i];
332           re = a->i[i+1];
333 
334           while (rs < re)
335             {
336               if (a->a[rs] != 0.0)
337                 {
338                   lusol->indc[nz] = i + 1;
339                   lusol->indr[nz] = a->j[rs] + 1;
340                   lusol->data[nz] = a->a[rs];
341                   nz++;
342                 }
343               rs++;
344             }
345         }
346 
347       /*******************************************************************/
348       /* Do the factorization.                                           */
349       /*******************************************************************/
350 
351       LU1FAC(&m, &n, &nz, &nnz,
352              lusol->luparm, lusol->parmlu, lusol->data,
353              lusol->indc, lusol->indr, lusol->ip, lusol->iq,
354              lusol->lenc, lusol->lenr, lusol->locc, lusol->locr,
355              lusol->iploc, lusol->iqloc, lusol->ipinv,
356              lusol->iqinv, lusol->mnsw, &status);
357 
358       switch(status)
359         {
360         case 0:		/* factored */
361           break;
362 
363         case 7:		/* insufficient memory */
364           break;
365 
366         case 1:
367         case -1:		/* singular */
368           SETERRQ(PETSC_ERR_LIB,"Singular matrix");
369 
370         case 3:
371         case 4:		/* error conditions */
372           SETERRQ(PETSC_ERR_LIB,"matrix error");
373 
374         default:		/* unknown condition */
375           SETERRQ(PETSC_ERR_LIB,"matrix unknown return code");
376         }
377 
378       factorizations++;
379     } while (status == 7);
380   (*F)->assembled = PETSC_TRUE;
381   PetscFunctionReturn(0);
382 }
383 
384 #undef __FUNCT__
385 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL"
386 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) {
387   /************************************************************************/
388   /* Input                                                                */
389   /*     A  - matrix to factor                                            */
390   /*     r  - row permutation (ignored)                                   */
391   /*     c  - column permutation (ignored)                                */
392   /*                                                                      */
393   /* Output                                                               */
394   /*     F  - matrix storing the factorization;                           */
395   /************************************************************************/
396   Mat       B;
397   Mat_LUSOL *lusol;
398   PetscErrorCode ierr;
399   int        i, m, n, nz, nnz;
400 
401   PetscFunctionBegin;
402 
403   /************************************************************************/
404   /* Check the arguments.                                                 */
405   /************************************************************************/
406 
407   ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr);
408   nz = ((Mat_SeqAIJ *)A->data)->nz;
409 
410   /************************************************************************/
411   /* Create the factorization.                                            */
412   /************************************************************************/
413 
414   ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,m,n,&B);CHKERRQ(ierr);
415   ierr = MatSetType(B,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 MatConvert_SeqAIJ_LUSOL(Mat A,const MatType type,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 (B != A) {
495     ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
496   }
497 
498   ierr                       = PetscNew(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   PetscLogInfo(0,"Using LUSOL for LU factorization and solves.");
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 MatCreate_LUSOL(Mat A)
557 {
558   PetscErrorCode ierr;
559 
560   PetscFunctionBegin;
561   /* Change type name before calling MatSetType to force proper construction of SeqAIJ and LUSOL types */
562   ierr = PetscObjectChangeTypeName((PetscObject)A,MATLUSOL);CHKERRQ(ierr);
563   ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr);
564   ierr = MatConvert_SeqAIJ_LUSOL(A,MATLUSOL,&A);CHKERRQ(ierr);
565   PetscFunctionReturn(0);
566 }
567 EXTERN_C_END
568