xref: /petsc/src/mat/impls/aij/seq/lusol/lusol.c (revision 69bb7ac94e38f481a514b28cb10ff7ece56113f9)
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     ierr = PetscFree(lusol->indc);CHKERRQ(ierr);
235   }
236 
237   ierr = MatConvert_LUSOL_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);
238   ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
239   PetscFunctionReturn(0);
240 }
241 
242 #undef __FUNCT__
243 #define __FUNCT__  "MatSolve_LUSOL"
244 PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x)
245 {
246   Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr;
247   double    *bb,*xx;
248   int       mode=5;
249   PetscErrorCode ierr;
250   int       i,m,n,nnz,status;
251 
252   PetscFunctionBegin;
253   ierr = VecGetArray(x, &xx);CHKERRQ(ierr);
254   ierr = VecGetArray(b, &bb);CHKERRQ(ierr);
255 
256   m = n = lusol->n;
257   nnz = lusol->nnz;
258 
259   for (i = 0; i < m; i++)
260     {
261       lusol->mnsv[i] = bb[i];
262     }
263 
264   LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz,
265          lusol->luparm, lusol->parmlu, lusol->data,
266          lusol->indc, lusol->indr, lusol->ip, lusol->iq,
267          lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status);
268 
269   if (status != 0)
270     {
271       SETERRQ(PETSC_ERR_ARG_SIZ,"solve failed");
272     }
273 
274   ierr = VecRestoreArray(x, &xx);CHKERRQ(ierr);
275   ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr);
276   PetscFunctionReturn(0);
277 }
278 
279 #undef __FUNCT__
280 #define __FUNCT__ "MatLUFactorNumeric_LUSOL"
281 PetscErrorCode MatLUFactorNumeric_LUSOL(Mat A,MatFactorInfo *info,Mat *F)
282 {
283   Mat_SeqAIJ     *a;
284   Mat_LUSOL      *lusol = (Mat_LUSOL*)(*F)->spptr;
285   PetscErrorCode ierr;
286   int            m, n, nz, nnz, status;
287   int            i, rs, re;
288   int            factorizations;
289 
290   PetscFunctionBegin;
291   ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);CHKERRQ(ierr);
292   a = (Mat_SeqAIJ *)A->data;
293 
294   if (m != lusol->n) {
295     SETERRQ(PETSC_ERR_ARG_SIZ,"factorization struct inconsistent");
296   }
297 
298   factorizations = 0;
299   do
300     {
301       /*******************************************************************/
302       /* Check the workspace allocation.                                 */
303       /*******************************************************************/
304 
305       nz = a->nz;
306       nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz));
307       nnz = PetscMax(nnz, 5*n);
308 
309       if (nnz < lusol->luparm[12]){
310         nnz = (int)(lusol->luroom * lusol->luparm[12]);
311       } else if ((factorizations > 0) && (lusol->luroom < 6)){
312         lusol->luroom += 0.1;
313       }
314 
315       nnz = PetscMax(nnz, (int)(lusol->luroom*(lusol->luparm[22] + lusol->luparm[23])));
316 
317       if (nnz > lusol->nnz){
318         ierr = PetscFree(lusol->indc);CHKERRQ(ierr);
319         ierr        = PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);CHKERRQ(ierr);
320         lusol->indr = lusol->indc + nnz;
321         lusol->data = (double *)(lusol->indr + nnz);
322         lusol->nnz  = nnz;
323       }
324 
325       /*******************************************************************/
326       /* Fill in the data for the problem.      (1-based Fortran style)  */
327       /*******************************************************************/
328 
329       nz = 0;
330       for (i = 0; i < n; i++)
331         {
332           rs = a->i[i];
333           re = a->i[i+1];
334 
335           while (rs < re)
336             {
337               if (a->a[rs] != 0.0)
338                 {
339                   lusol->indc[nz] = i + 1;
340                   lusol->indr[nz] = a->j[rs] + 1;
341                   lusol->data[nz] = a->a[rs];
342                   nz++;
343                 }
344               rs++;
345             }
346         }
347 
348       /*******************************************************************/
349       /* Do the factorization.                                           */
350       /*******************************************************************/
351 
352       LU1FAC(&m, &n, &nz, &nnz,
353              lusol->luparm, lusol->parmlu, lusol->data,
354              lusol->indc, lusol->indr, lusol->ip, lusol->iq,
355              lusol->lenc, lusol->lenr, lusol->locc, lusol->locr,
356              lusol->iploc, lusol->iqloc, lusol->ipinv,
357              lusol->iqinv, lusol->mnsw, &status);
358 
359       switch(status)
360         {
361         case 0:		/* factored */
362           break;
363 
364         case 7:		/* insufficient memory */
365           break;
366 
367         case 1:
368         case -1:		/* singular */
369           SETERRQ(PETSC_ERR_LIB,"Singular matrix");
370 
371         case 3:
372         case 4:		/* error conditions */
373           SETERRQ(PETSC_ERR_LIB,"matrix error");
374 
375         default:		/* unknown condition */
376           SETERRQ(PETSC_ERR_LIB,"matrix unknown return code");
377         }
378 
379       factorizations++;
380     } while (status == 7);
381   (*F)->assembled = PETSC_TRUE;
382   PetscFunctionReturn(0);
383 }
384 
385 #undef __FUNCT__
386 #define __FUNCT__ "MatLUFactorSymbolic_LUSOL"
387 PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) {
388   /************************************************************************/
389   /* Input                                                                */
390   /*     A  - matrix to factor                                            */
391   /*     r  - row permutation (ignored)                                   */
392   /*     c  - column permutation (ignored)                                */
393   /*                                                                      */
394   /* Output                                                               */
395   /*     F  - matrix storing the factorization;                           */
396   /************************************************************************/
397   Mat       B;
398   Mat_LUSOL *lusol;
399   PetscErrorCode ierr;
400   int        i, m, n, nz, nnz;
401 
402   PetscFunctionBegin;
403 
404   /************************************************************************/
405   /* Check the arguments.                                                 */
406   /************************************************************************/
407 
408   ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr);
409   nz = ((Mat_SeqAIJ *)A->data)->nz;
410 
411   /************************************************************************/
412   /* Create the factorization.                                            */
413   /************************************************************************/
414 
415   ierr = MatCreate(A->comm,&B);CHKERRQ(ierr);
416   ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);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 = PetscInfo(0,"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