IFPACK Development
Loading...
Searching...
No Matches
Ifpack_IC_Utils.cpp
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
63void 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
70void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
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
80static 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
135static 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
205static 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
240static void update_lists_newcol(
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
275void crout_ict(
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,
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}