Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_IC_Utils.cpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "Ifpack_ConfigDefs.h"
44 #include "Ifpack_IC_Utils.h"
45 
46 #define SYMSTR 1 /* structurally symmetric version */
47 #include <stdio.h>
48 
49 #define MIN(a,b) ((a)<=(b) ? (a) : (b))
50 #define MAX(a,b) ((a)>=(b) ? (a) : (b))
51 #define ABS(a) ((a)>=0 ? (a) : -(a))
52 
53 #define SHORTCUT(p, a, ja, ia) \
54  (a) = (p)->val; \
55  (ja) = (p)->col; \
56  (ia) = (p)->ptr;
57 
58 #define MATNULL(p) \
59  (p).val = NULL; \
60  (p).col = NULL; \
61  (p).ptr = NULL;
62 
63 void Ifpack_AIJMatrix_alloc(Ifpack_AIJMatrix *a, int n, int nnz)
64 {
65  a->val = new double[nnz];
66  a->col = new int[nnz];
67  a->ptr = new int[n+1];
68 }
69 
71 {
72  delete [] (a->val);
73  delete [] (a->col);
74  delete [] (a->ptr);
75  a->val = 0;
76  a->col = 0;
77  a->ptr = 0;
78 }
79 
80 static void qsplit(double *a, int *ind, int n, int ncut)
81 {
82  double tmp, abskey;
83  int itmp, first, last, mid;
84  int j;
85 
86  ncut--;
87  first = 0;
88  last = n-1;
89  if (ncut < first || ncut > last)
90  return;
91 
92  /* outer loop while mid != ncut */
93  while (1)
94  {
95  mid = first;
96  abskey = ABS(a[mid]);
97  for (j=first+1; j<=last; j++)
98  {
99  if (ABS(a[j]) > abskey)
100  {
101  mid = mid+1;
102  /* interchange */
103  tmp = a[mid];
104  itmp = ind[mid];
105  a[mid] = a[j];
106  ind[mid] = ind[j];
107  a[j] = tmp;
108  ind[j] = itmp;
109  }
110  }
111 
112  /* interchange */
113  tmp = a[mid];
114  a[mid] = a[first];
115  a[first] = tmp;
116  itmp = ind[mid];
117  ind[mid] = ind[first];
118  ind[first] = itmp;
119 
120  /* test for while loop */
121  if (mid == ncut)
122  return;
123 
124  if (mid > ncut)
125  last = mid-1;
126  else
127  first = mid+1;
128  }
129 }
130 
131 /* update column k using previous columns */
132 /* assumes that column of A resides in the work vector */
133 /* this function can also be used to update rows */
134 
135 static void update_column(
136  int k,
137  const int *ia, const int *ja, const double *a,
138  const int *ifirst,
139  const int *ifirst2,
140  const int *list2,
141  const double *multipliers, /* the val array of the other factor */
142  const double *d, /* diagonal of factorization */
143  int *marker,
144  double *ta,
145  int *itcol,
146  int *ptalen)
147 {
148  int j, i, isj, iej, row;
149  int talen, pos;
150  double mult;
151 
152  talen = *ptalen;
153 
154  j = list2[k];
155  while (j != -1)
156  {
157  /* update column k using column j */
158 
159  isj = ifirst[j];
160 
161  /* skip first nonzero in column, since it does not contribute */
162  /* if symmetric structure */
163  /* isj++; */
164 
165  /* now do the actual update */
166  if (isj != -1)
167  {
168  /* multiplier */
169  mult = multipliers[ifirst2[j]] * d[j];
170 
171  /* section of column used for update */
172  iej = ia[j+1]-1;
173 
174  for (i=isj; i<=iej; i++)
175  {
176  row = ja[i];
177 
178  /* if nonsymmetric structure */
179  if (row <= k)
180  continue;
181 
182  if ((pos = marker[row]) != -1)
183  {
184  /* already in pattern of working row */
185  ta[pos] -= mult*a[i];
186  }
187  else
188  {
189  /* not yet in pattern of working row */
190  itcol[talen] = row;
191  ta[talen] = -mult*a[i];
192  marker[row] = talen++;
193  }
194  }
195  }
196 
197  j = list2[j];
198  }
199 
200  *ptalen = talen;
201 }
202 
203 /* update ifirst and list */
204 
205 static void update_lists(
206  int k,
207  const int *ia,
208  const int *ja,
209  int *ifirst,
210  int *list)
211 {
212  int j, isj, iej, iptr;
213 
214  j = list[k];
215  while (j != -1)
216  {
217  isj = ifirst[j];
218  iej = ia[j+1]-1;
219 
220  isj++;
221 
222  if (isj != 0 && isj <= iej) /* nonsymmetric structure */
223  {
224  /* increment ifirst for column j */
225  ifirst[j] = isj;
226 
227  /* add j to head of list for list[ja[isj]] */
228  iptr = j;
229  j = list[j];
230  list[iptr] = list[ja[isj]];
231  list[ja[isj]] = iptr;
232  }
233  else
234  {
235  j = list[j];
236  }
237  }
238 }
239 
241  int k,
242  int isk,
243  int iptr,
244  int *ifirst,
245  int *list)
246 {
247  ifirst[k] = isk;
248  list[k] = list[iptr];
249  list[iptr] = k;
250 }
251 
252 /*
253  * crout_ict - Crout version of ICT - Incomplete Cholesky by Threshold
254  *
255  * The incomplete factorization L D L^T is computed,
256  * where L is unit triangular, and D is diagonal
257  *
258  * INPUTS
259  * n = number of rows or columns of the matrices
260  * AL = unit lower triangular part of A stored by columns
261  * the unit diagonal is implied (not stored)
262  * Adiag = diagonal part of A
263  * droptol = drop tolerance
264  * lfil = max nonzeros per col in L factor or per row in U factor
265  * TODO: lfil should rather be a fill ratio applied to each column
266  *
267  * OUTPUTS
268  * L = lower triangular factor stored by columns
269  * pdiag = diagonal factor stored in an array
270  *
271  * NOTE: calling function must free the memory allocated by this
272  * function, i.e., L, pdiag.
273  */
274 
276  int n,
277 #ifdef IFPACK
278  void * A,
279  int maxentries;
280  int (*getcol)( void * A, int col, int ** nentries, double * val, int * ind),
281  int (*getdiag)( void *A, double * diag),
282 #else
283  const Ifpack_AIJMatrix *AL,
284  const double *Adiag,
285 #endif
286  double droptol,
287  int lfil,
288  Ifpack_AIJMatrix *L,
289  double **pdiag)
290 {
291  int k, j, i, index;
292  int count_l;
293  double norm_l;
294 
295  /* work arrays; work_l is list of nonzero values */
296  double *work_l = new double[n];
297  int *ind_l = new int[n];
298  int len_l;
299 
300  /* list and ifirst data structures */
301  int *list_l = new int[n];
302  int *ifirst_l = new int[n];
303 
304  /* aliases */
305  int *marker_l = ifirst_l;
306 
307  /* matrix arrays */
308  double *al; int *jal, *ial;
309  double *l; int *jl, *il;
310 
311  int kl = 0;
312 
313  double *diag = new double[n];
314  *pdiag = diag;
315 
316  Ifpack_AIJMatrix_alloc(L, n, lfil*n);
317 
318  SHORTCUT(AL, al, jal, ial);
319  SHORTCUT(L, l, jl, il);
320 
321  /* initialize work arrays */
322  for (i=0; i<n; i++)
323  {
324  list_l[i] = -1;
325  ifirst_l[i] = -1;
326  marker_l[i] = -1;
327  }
328 
329  /* copy the diagonal */
330 #ifdef IFPACK
331  getdiag(A, diag);
332 #else
333  for (i=0; i<n; i++)
334  diag[i] = Adiag[i];
335 #endif
336 
337  /* start off the matrices right */
338  il[0] = 0;
339 
340  /*
341  * Main loop over columns and rows
342  */
343 
344  for (k=0; k<n; k++)
345  {
346  /*
347  * lower triangular factor update by columns
348  * (need ifirst for L and lists for U)
349  */
350 
351  /* copy column of A into work vector */
352  norm_l = 0.;
353 #ifdef IFPACK
354  getcol(A, k, len_l, work_l, ind_l);
355  for (j=0; j<len_l; j++)
356  {
357  norm_l += ABS(work_l[j]);
358  marker_l[ind_l[j]] = j;
359  }
360 #else
361  len_l = 0;
362  for (j=ial[k]; j<ial[k+1]; j++)
363  {
364  index = jal[j];
365  work_l[len_l] = al[j];
366  norm_l += ABS(al[j]);
367  ind_l[len_l] = index;
368  marker_l[index] = len_l++; /* points into work array */
369  }
370 #endif
371  norm_l = (len_l == 0) ? 0.0 : norm_l/((double) len_l);
372 
373  /* update and scale */
374 
375  update_column(k, il, jl, l, ifirst_l, ifirst_l, list_l, l,
376  diag, marker_l, work_l, ind_l, &len_l);
377 
378  for (j=0; j<len_l; j++)
379  work_l[j] /= diag[k];
380 
381  i = 0;
382  for (j=0; j<len_l; j++)
383  {
384  if (ABS(work_l[j]) < droptol * norm_l)
385  {
386  /* zero out marker array for this */
387  marker_l[ind_l[j]] = -1;
388  }
389  else
390  {
391  work_l[i] = work_l[j];
392  ind_l[i] = ind_l[j];
393  i++;
394  }
395  }
396  len_l = i;
397 
398  /*
399  * dropping: for each work vector, perform qsplit, and then
400  * sort each part by increasing index; then copy each sorted
401  * part into the factors
402  */
403 
404  // TODO: keep different #nonzeros per column depending on original matrix
405  count_l = MIN(len_l, lfil);
406  qsplit(work_l, ind_l, len_l, count_l);
407  ifpack_multilist_sort(ind_l, work_l, count_l);
408 
409  for (j=0; j<count_l; j++)
410  {
411  l[kl] = work_l[j];
412  jl[kl++] = ind_l[j];
413  }
414  il[k+1] = kl;
415 
416  /*
417  * update lists
418  */
419 
420  update_lists(k, il, jl, ifirst_l, list_l);
421 
422  if (kl - il[k] > SYMSTR)
423  update_lists_newcol(k, il[k], jl[il[k]], ifirst_l, list_l);
424 
425  /* zero out the marker arrays */
426  for (j=0; j<len_l; j++)
427  marker_l[ind_l[j]] = -1;
428 
429  /*
430  * update the diagonal (after dropping)
431  */
432 
433  for (j=0; j<count_l; j++)
434  {
435  index = ind_l[j];
436  diag[index] -= work_l[j] * work_l[j] * diag[k];
437  }
438  }
439 
440  delete [] work_l;
441  delete [] ind_l;
442  delete [] list_l;
443  delete [] ifirst_l;
444 }
#define SHORTCUT(p, a, ja, ia)
static void update_lists(int k, const int *ia, const int *ja, int *ifirst, int *list)
void Ifpack_AIJMatrix_alloc(Ifpack_AIJMatrix *a, int n, int nnz)
void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
static void update_lists_newcol(int k, int isk, int iptr, int *ifirst, int *list)
static void qsplit(double *a, int *ind, int n, int ncut)
#define ABS(a)
#define SYMSTR
void ifpack_multilist_sort(int *const pbase, double *const daux, int total_elems)
static void update_column(int k, const int *ia, const int *ja, const double *a, const int *ifirst, const int *ifirst2, const int *list2, const double *multipliers, const double *d, int *marker, double *ta, int *itcol, int *ptalen)
void crout_ict(int n, const Ifpack_AIJMatrix *AL, const double *Adiag, double droptol, int lfil, Ifpack_AIJMatrix *L, double **pdiag)
#define MIN(a, b)