superlu_complex.c
Go to the documentation of this file.
1
2 /* Wrapper for SuperLU solver. Based on fortran wrapper supplied
3  * with SuperLU version 3.0. Given that, it seems appropriate
4  * to retain this copyright notice:
5  *
6  * -- SuperLU routine (version 3.0) --
7  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
8  * and Lawrence Berkeley National Lab.
9  * October 15, 2003
10  *
11  */
12
13 #include "../../external_src/oomph_superlu_4.3/slu_zdefs.h"
14 #include "math.h"
15
16
17 /* ================================================= */
18 /* Pointer to the LU factors*/
19 /* ================================================= */
20 typedef void* fptr;
21
22
23 /* ================================================= */
24 /* Struct for the lu factors */
25 /* ================================================= */
26 typedef struct {
27  SuperMatrix *L;
28  SuperMatrix *U;
29  int *perm_c;
30  int *perm_r;
31 } factors_t;
32
33
34
35 /* =========================================================================
36  * Wrapper to superlu solver:
37  *
38  * op_flag = int specifies the operation:
39  * 1, performs LU decomposition for the first time
40  * 2, performs triangular solve
41  * 3, free all the storage in the end
42  * n = dimension of matrix
43  * nnz = # of nonzero entries
44  * nrhs = # of RHSs
45  * values = double array containing the nonzero entries
46  * rowind = int array containing the row indices of the entries
47  * colptr = int array containing the column start
48  * b = double array containing the rhs vector (gets overwritten
49  * with solution)
50  * ldb = leading dimension of b
51  * transpose = 0/1 if matrix is transposed/not transposed
52  * doc = 0/1 for full doc/no full doc
53  * info = info flag from superlu
54  * f_factors = pointer to LU factors. (If op_flag == 1, it is an output
55  * and contains the pointer pointing to the structure of
56  * the factored matrices. Otherwise, it it an input.
57  * Returns the SIGN of the determinant of the matrix
58  * =========================================================================
59  */
60 int superlu_complex(int *op_flag, int *n, int *nnz, int *nrhs,
61  doublecomplex *values, int *rowind, int *colptr,
62  doublecomplex *b, int *ldb, int *transpose, int *doc,
63  fptr *f_factors, int *info)
64
65 {
66
67  SuperMatrix A, AC, B;
68  SuperMatrix *L, *U;
69  int *perm_r; /* row permutations from partial pivoting */
70  int *perm_c; /* column permutation vector */
71  int *etree; /* column elimination tree */
72  SCformat *Lstore;
73  NCformat *Ustore;
74  int i, j, panel_size, permc_spec, relax;
75  trans_t trans;
76  //double drop_tol = 0.0; //No longer needed SuperLU 4.3
77  mem_usage_t mem_usage;
78  superlu_options_t options;
79  SuperLUStat_t stat;
80  factors_t *LUfactors;
81
82  doublecomplex *Lval;
83  doublecomplex *diagU, *dblock;
84  int_t fsupc, nsupr, nsupc, luptr;
85  int_t i2, k2, nsupers;
86  int signature=1;
87  int sign = 0;
88
89 /* Do we need to transpose? */
90  if (*transpose==0)
91  {
92  trans = NOTRANS;
93  }
94  else
95  {
96  trans = TRANS;
97  }
98
99  if ( *op_flag == 1 ) { /* LU decomposition */
100
101  /* Set the default input options. */
102  set_default_options(&options);
103
104  /* Initialize the statistics variables. */
105  StatInit(&stat);
106
107  zCreate_CompCol_Matrix(&A, *n, *n, *nnz, values, rowind, colptr,
108  SLU_NC, SLU_Z, SLU_GE);
109  L = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
110  U = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
111  if ( !(perm_r = intMalloc(*n)) ) ABORT("Malloc fails for perm_r[].");
112  if ( !(perm_c = intMalloc(*n)) ) ABORT("Malloc fails for perm_c[].");
113  if ( !(etree = intMalloc(*n)) ) ABORT("Malloc fails for etree[].");
114
115  /*
116  * Get column permutation vector perm_c[], according to permc_spec:
117  * permc_spec = 0: natural ordering
118  * permc_spec = 1: minimum degree on structure of A'*A
119  * permc_spec = 2: minimum degree on structure of A'+A
120  * permc_spec = 3: approximate minimum degree for unsymmetric matrices
121  */
122  permc_spec = options.ColPerm;
123  get_perm_c(permc_spec, &A, perm_c);
124
125  sp_preorder(&options, &A, perm_c, etree, &AC);
126
127  panel_size = sp_ienv(1);
128  relax = sp_ienv(2);
129
130  zgstrf(&options, &AC, /*drop_tol,*/ relax, panel_size,
131  etree, NULL, 0, perm_c, perm_r, L, U, &stat, info);
132
133  if ( *info == 0 ) {
134  Lstore = (SCformat *) L->Store;
135  Ustore = (NCformat *) U->Store;
136  if (*doc!=0)
137  {
138  printf(" No of nonzeros in factor L = %d\n", Lstore->nnz);
139  printf(" No of nonzeros in factor U = %d\n", Ustore->nnz);
140  printf(" No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz);
141  zQuerySpace(L, U, &mem_usage);
142  printf(" L\\U MB %.3f\ttotal MB needed %.3f\n",
143  //\texpansions %d\n",
144  mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
145  //mem_usage.expansions);
146  }
147  } else {
148  printf("dgstrf() error returns INFO= %d\n", *info);
149  if ( *info <= *n ) { /* factorization completes */
150  zQuerySpace(L, U, &mem_usage);
151  printf(" L\\U MB %.3f\ttotal MB needed %.3f\n",
152  //\texpansions %d\n\n",
153  mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
154  //mem_usage.expansions);
155  }
156  }
157
158  /* Save the LU factors in the factors handle */
159  LUfactors = (factors_t*) SUPERLU_MALLOC(sizeof(factors_t));
160  LUfactors->L = L;
161  LUfactors->U = U;
162  LUfactors->perm_c = perm_c;
163  LUfactors->perm_r = perm_r;
164  *f_factors = (fptr) LUfactors;
165
166  //Work out and print the sign of the determinant
167  //This code is hacked from supraLU by Alex Pletzer
168  //and Doug McCune (NTCC) (http://w3.pppl.gov/ntcc/SupraLu/)
169  Lstore = L->Store;
170  Lval = Lstore->nzval;
171  nsupers = Lstore->nsuper + 1;
172
173  //Get the diagonal entries of the U matrix
174  //Allocate store for the entries
175  if ( !(diagU = SUPERLU_MALLOC( *n * sizeof(SuperMatrix) )) )
176  ABORT("Malloc fails for diagU[].");
177  //Loop over the number of super diagonal terms(?)
178  for(k2=0; k2< nsupers; k2++)
179  {
180  fsupc = L_FST_SUPC(k2);
181  nsupc = L_FST_SUPC(k2+1) - fsupc;
182  nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);
183  luptr = L_NZ_START(fsupc);
184
185  dblock = &diagU[fsupc];
186  for(i2 = 0; i2 < nsupc; i2++)
187  {
188  dblock[i2] = Lval[luptr];
189  luptr += nsupr + 1;
190  }
191  }
192
193  //Now multiply all the diagonal terms together to get the determinant
194  //Note that we need to use the mantissa, exponent formulation to
195  //avoid underflow errors
196  double determinant_mantissa=1.0;
197  int determinant_exponent = 0, iexp;
198  for(i=0; i<*n; i++)
199  {
200  determinant_mantissa *= frexp(diagU[i].r, &iexp);
201  determinant_exponent += iexp;
202  /* normalise*/
203  determinant_mantissa = frexp(determinant_mantissa,&iexp);
204  determinant_exponent += iexp;
205
206  /*Now worry about the permutations
207  (this is done in a stupid, but not too inefficient way)*/
208  for(j=i;j<*n;j++)
209  {
210  if(perm_r[j] < perm_r[i]) {signature *= -1;}
211  if(perm_c[j] < perm_c[i]) {signature *= -1;}
212  }
213  }
214
215  //Find the sign of the determinant
216  if(determinant_mantissa > 0.0) {sign = 1;}
217  if(determinant_mantissa < 0.0) {sign = -1;}
218
219  //Multiply the sign by the signature
220  sign *= signature;
221
222  /* Free un-wanted storage */
223  SUPERLU_FREE(diagU);
224  SUPERLU_FREE(etree);
225  Destroy_SuperMatrix_Store(&A);
226  Destroy_CompCol_Permuted(&AC);
227  StatFree(&stat);
228
229  //Return the sign of the determinant
230  return sign;
231
232  } else if ( *op_flag == 2 ) { /* Triangular solve */
233  /* Initialize the statistics variables. */
234  StatInit(&stat);
235
236  /* Extract the LU factors in the factors handle */
237  LUfactors = (factors_t*) *f_factors;
238  L = LUfactors->L;
239  U = LUfactors->U;
240  perm_c = LUfactors->perm_c;
241  perm_r = LUfactors->perm_r;
242
243  zCreate_Dense_Matrix(&B, *n, *nrhs, b, *ldb, SLU_DN, SLU_Z, SLU_GE);
244
245  /* Solve the system A*X=B, overwriting B with X. */
246  zgstrs (trans, L, U, perm_c, perm_r, &B, &stat, info);
247
248  Destroy_SuperMatrix_Store(&B);
249  StatFree(&stat);
250
251  //Return zero
252  return 0;
253
254  } else if ( *op_flag == 3 ) { /* Free storage */
255  /* Free the LU factors in the factors handle */
256  LUfactors = (factors_t*) *f_factors;
257  SUPERLU_FREE (LUfactors->perm_r);
258  SUPERLU_FREE (LUfactors->perm_c);
259  Destroy_SuperNode_Matrix(LUfactors->L);
260  Destroy_CompCol_Matrix(LUfactors->U);
261  SUPERLU_FREE (LUfactors->L);
262  SUPERLU_FREE (LUfactors->U);
263  SUPERLU_FREE (LUfactors);
264  return 0;
265  } else {
266  fprintf(stderr,"Invalid op_flag=%d passed to c_cpp_dgssv()\n",*op_flag);
267  exit(-1);
268  }
269 }
270
271
int * perm_c
Definition: superlu.c:32
cstr elem_len * i
Definition: cfortran.h:607
int superlu_complex(int *op_flag, int *n, int *nnz, int *nrhs, doublecomplex *values, int *rowind, int *colptr, doublecomplex *b, int *ldb, int *transpose, int *doc, fptr *f_factors, int *info)
int * perm_r
Definition: superlu.c:33
void * fptr
void * fptr
Definition: superlu.c:23
SuperMatrix * U
Definition: superlu.c:31
SuperMatrix * L
Definition: superlu.c:30