00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_IC_Utils.h"
00003
00004 #define SYMSTR 1
00005 #include <stdio.h>
00006
00007 #define MIN(a,b) ((a)<=(b) ? (a) : (b))
00008 #define MAX(a,b) ((a)>=(b) ? (a) : (b))
00009 #define ABS(a) ((a)>=0 ? (a) : -(a))
00010
00011 #define SHORTCUT(p, a, ja, ia) \
00012 (a) = (p)->val; \
00013 (ja) = (p)->col; \
00014 (ia) = (p)->ptr;
00015
00016 #define MATNULL(p) \
00017 (p).val = NULL; \
00018 (p).col = NULL; \
00019 (p).ptr = NULL;
00020
00021 void Ifpack_AIJMatrix_alloc(Ifpack_AIJMatrix *a, int n, int nnz)
00022 {
00023 a->val = new double[nnz];
00024 a->col = new int[nnz];
00025 a->ptr = new int[n+1];
00026 }
00027
00028 void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
00029 {
00030 delete [] (a->val);
00031 delete [] (a->col);
00032 delete [] (a->ptr);
00033 a->val = 0;
00034 a->col = 0;
00035 a->ptr = 0;
00036 }
00037
00038 static void qsplit(double *a, int *ind, int n, int ncut)
00039 {
00040 double tmp, abskey;
00041 int itmp, first, last, mid;
00042 int j;
00043
00044 ncut--;
00045 first = 0;
00046 last = n-1;
00047 if (ncut < first || ncut > last)
00048 return;
00049
00050
00051 while (1)
00052 {
00053 mid = first;
00054 abskey = ABS(a[mid]);
00055 for (j=first+1; j<=last; j++)
00056 {
00057 if (ABS(a[j]) > abskey)
00058 {
00059 mid = mid+1;
00060
00061 tmp = a[mid];
00062 itmp = ind[mid];
00063 a[mid] = a[j];
00064 ind[mid] = ind[j];
00065 a[j] = tmp;
00066 ind[j] = itmp;
00067 }
00068 }
00069
00070
00071 tmp = a[mid];
00072 a[mid] = a[first];
00073 a[first] = tmp;
00074 itmp = ind[mid];
00075 ind[mid] = ind[first];
00076 ind[first] = itmp;
00077
00078
00079 if (mid == ncut)
00080 return;
00081
00082 if (mid > ncut)
00083 last = mid-1;
00084 else
00085 first = mid+1;
00086 }
00087 }
00088
00089
00090
00091
00092
00093 static void update_column(
00094 int k,
00095 const int *ia, const int *ja, const double *a,
00096 const int *ifirst,
00097 const int *ifirst2,
00098 const int *list2,
00099 const double *multipliers,
00100 const double *d,
00101 int *marker,
00102 double *ta,
00103 int *itcol,
00104 int *ptalen)
00105 {
00106 int j, i, isj, iej, row;
00107 int talen, pos;
00108 double mult;
00109
00110 talen = *ptalen;
00111
00112 j = list2[k];
00113 while (j != -1)
00114 {
00115
00116
00117 isj = ifirst[j];
00118
00119
00120
00121
00122
00123
00124 if (isj != -1)
00125 {
00126
00127 mult = multipliers[ifirst2[j]] * d[j];
00128
00129
00130 iej = ia[j+1]-1;
00131
00132 for (i=isj; i<=iej; i++)
00133 {
00134 row = ja[i];
00135
00136
00137 if (row <= k)
00138 continue;
00139
00140 if ((pos = marker[row]) != -1)
00141 {
00142
00143 ta[pos] -= mult*a[i];
00144 }
00145 else
00146 {
00147
00148 itcol[talen] = row;
00149 ta[talen] = -mult*a[i];
00150 marker[row] = talen++;
00151 }
00152 }
00153 }
00154
00155 j = list2[j];
00156 }
00157
00158 *ptalen = talen;
00159 }
00160
00161
00162
00163 static void update_lists(
00164 int k,
00165 const int *ia,
00166 const int *ja,
00167 int *ifirst,
00168 int *list)
00169 {
00170 int j, isj, iej, iptr;
00171
00172 j = list[k];
00173 while (j != -1)
00174 {
00175 isj = ifirst[j];
00176 iej = ia[j+1]-1;
00177
00178 isj++;
00179
00180 if (isj != 0 && isj <= iej)
00181 {
00182
00183 ifirst[j] = isj;
00184
00185
00186 iptr = j;
00187 j = list[j];
00188 list[iptr] = list[ja[isj]];
00189 list[ja[isj]] = iptr;
00190 }
00191 else
00192 {
00193 j = list[j];
00194 }
00195 }
00196 }
00197
00198 static void update_lists_newcol(
00199 int k,
00200 int isk,
00201 int iptr,
00202 int *ifirst,
00203 int *list)
00204 {
00205 ifirst[k] = isk;
00206 list[k] = list[iptr];
00207 list[iptr] = k;
00208 }
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232 void crout_ict(
00233 int n,
00234 #ifdef IFPACK
00235 void * A,
00236 int maxentries;
00237 int (*getcol)( void * A, int col, int ** nentries, double * val, int * ind),
00238 int (*getdiag)( void *A, double * diag),
00239 #else
00240 const Ifpack_AIJMatrix *AL,
00241 const double *Adiag,
00242 #endif
00243 double droptol,
00244 int lfil,
00245 Ifpack_AIJMatrix *L,
00246 double **pdiag)
00247 {
00248 int k, j, i, index;
00249 int count_l;
00250 double norm_l;
00251
00252
00253 double *work_l = new double[n];
00254 int *ind_l = new int[n];
00255 int len_l;
00256
00257
00258 int *list_l = new int[n];
00259 int *ifirst_l = new int[n];
00260
00261
00262 int *marker_l = ifirst_l;
00263
00264
00265 double *al; int *jal, *ial;
00266 double *l; int *jl, *il;
00267
00268 int kl = 0;
00269
00270 double *diag = new double[n];
00271 *pdiag = diag;
00272
00273 Ifpack_AIJMatrix_alloc(L, n, lfil*n);
00274
00275 SHORTCUT(AL, al, jal, ial);
00276 SHORTCUT(L, l, jl, il);
00277
00278
00279 for (i=0; i<n; i++)
00280 {
00281 list_l[i] = -1;
00282 ifirst_l[i] = -1;
00283 marker_l[i] = -1;
00284 }
00285
00286
00287 #ifdef IFPACK
00288 getdiag(A, diag);
00289 #else
00290 for (i=0; i<n; i++)
00291 diag[i] = Adiag[i];
00292 #endif
00293
00294
00295 il[0] = 0;
00296
00297
00298
00299
00300
00301 for (k=0; k<n; k++)
00302 {
00303
00304
00305
00306
00307
00308
00309 norm_l = 0.;
00310 #ifdef IFPACK
00311 getcol(A, k, len_l, work_l, ind_l);
00312 for (j=0; j<len_l; j++)
00313 {
00314 norm_l += ABS(work_l[j]);
00315 marker_l[ind_l[j]] = j;
00316 }
00317 #else
00318 len_l = 0;
00319 for (j=ial[k]; j<ial[k+1]; j++)
00320 {
00321 index = jal[j];
00322 work_l[len_l] = al[j];
00323 norm_l += ABS(al[j]);
00324 ind_l[len_l] = index;
00325 marker_l[index] = len_l++;
00326 }
00327 #endif
00328 norm_l = (len_l == 0) ? 0.0 : norm_l/((double) len_l);
00329
00330
00331
00332 update_column(k, il, jl, l, ifirst_l, ifirst_l, list_l, l,
00333 diag, marker_l, work_l, ind_l, &len_l);
00334
00335 for (j=0; j<len_l; j++)
00336 work_l[j] /= diag[k];
00337
00338 i = 0;
00339 for (j=0; j<len_l; j++)
00340 {
00341 if (ABS(work_l[j]) < droptol * norm_l)
00342 {
00343
00344 marker_l[ind_l[j]] = -1;
00345 }
00346 else
00347 {
00348 work_l[i] = work_l[j];
00349 ind_l[i] = ind_l[j];
00350 i++;
00351 }
00352 }
00353 len_l = i;
00354
00355
00356
00357
00358
00359
00360
00361 count_l = MIN(len_l, lfil);
00362 qsplit(work_l, ind_l, len_l, count_l);
00363 quicksort(ind_l, work_l, count_l);
00364
00365 for (j=0; j<count_l; j++)
00366 {
00367 l[kl] = work_l[j];
00368 jl[kl++] = ind_l[j];
00369 }
00370 il[k+1] = kl;
00371
00372
00373
00374
00375
00376 update_lists(k, il, jl, ifirst_l, list_l);
00377
00378 if (kl - il[k] > SYMSTR)
00379 update_lists_newcol(k, il[k], jl[il[k]], ifirst_l, list_l);
00380
00381
00382 for (j=0; j<len_l; j++)
00383 marker_l[ind_l[j]] = -1;
00384
00385
00386
00387
00388
00389 for (j=0; j<count_l; j++)
00390 {
00391 index = ind_l[j];
00392 diag[index] -= work_l[j] * work_l[j] * diag[k];
00393 }
00394 }
00395
00396 delete [] work_l;
00397 delete [] ind_l;
00398 delete [] list_l;
00399 delete [] ifirst_l;
00400 }