| 1 |
|
| 2 |
/******************************************************* |
| 3 |
* |
| 4 |
* Copyright (c) 2003-2008 by University of Queensland |
| 5 |
* Earth Systems Science Computational Center (ESSCC) |
| 6 |
* http://www.uq.edu.au/esscc |
| 7 |
* |
| 8 |
* Primary Business: Queensland, Australia |
| 9 |
* Licensed under the Open Software License version 3.0 |
| 10 |
* http://www.opensource.org/licenses/osl-3.0.php |
| 11 |
* |
| 12 |
*******************************************************/ |
| 13 |
|
| 14 |
|
| 15 |
/**************************************************************/ |
| 16 |
|
| 17 |
/* Paso: AMG preconditioner with reordering */ |
| 18 |
|
| 19 |
/**************************************************************/ |
| 20 |
|
| 21 |
/* Author: artak@access.edu.au */ |
| 22 |
|
| 23 |
/**************************************************************/ |
| 24 |
|
| 25 |
#include "Paso.h" |
| 26 |
#include "Solver.h" |
| 27 |
#include "PasoUtil.h" |
| 28 |
#include "Pattern_coupling.h" |
| 29 |
|
| 30 |
/**************************************************************/ |
| 31 |
|
| 32 |
/* free all memory used by AMG */ |
| 33 |
|
| 34 |
void Paso_Solver_AMG_free(Paso_Solver_AMG * in) { |
| 35 |
if (in!=NULL) { |
| 36 |
Paso_Solver_AMG_free(in->AMG_of_Schur); |
| 37 |
MEMFREE(in->inv_A_FF); |
| 38 |
MEMFREE(in->A_FF_pivot); |
| 39 |
Paso_SparseMatrix_free(in->A_FC); |
| 40 |
Paso_SparseMatrix_free(in->A_CF); |
| 41 |
MEMFREE(in->rows_in_F); |
| 42 |
MEMFREE(in->rows_in_C); |
| 43 |
MEMFREE(in->mask_F); |
| 44 |
MEMFREE(in->mask_C); |
| 45 |
MEMFREE(in->x_F); |
| 46 |
MEMFREE(in->b_F); |
| 47 |
MEMFREE(in->x_C); |
| 48 |
MEMFREE(in->b_C); |
| 49 |
MEMFREE(in); |
| 50 |
} |
| 51 |
} |
| 52 |
|
| 53 |
/**************************************************************/ |
| 54 |
|
| 55 |
/* constructs the block-block factorization of |
| 56 |
|
| 57 |
[ A_FF A_FC ] |
| 58 |
A_p= |
| 59 |
[ A_CF A_FF ] |
| 60 |
|
| 61 |
to |
| 62 |
|
| 63 |
[ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ] |
| 64 |
[ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] |
| 65 |
|
| 66 |
|
| 67 |
where S=A_FF-ACF*invA_FF*A_FC within the shape of S |
| 68 |
|
| 69 |
then AMG is applied to S again until S becomes empty |
| 70 |
|
| 71 |
*/ |
| 72 |
Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,bool_t verbose,dim_t level) { |
| 73 |
Paso_Solver_AMG* out=NULL; |
| 74 |
dim_t n=A_p->numRows; |
| 75 |
dim_t n_block=A_p->row_block_size; |
| 76 |
index_t* mis_marker=NULL; |
| 77 |
index_t* counter=NULL; |
| 78 |
index_t iPtr,*index, *where_p, iPtr_s; |
| 79 |
dim_t i,k,j,j0; |
| 80 |
Paso_SparseMatrix * schur=NULL; |
| 81 |
Paso_SparseMatrix * schur_withFillIn=NULL; |
| 82 |
Paso_Pattern* mult=NULL; |
| 83 |
double time0,time1,time2,S; |
| 84 |
|
| 85 |
/* identify independend set of rows/columns */ |
| 86 |
mis_marker=TMPMEMALLOC(n,index_t); |
| 87 |
counter=TMPMEMALLOC(n,index_t); |
| 88 |
out=MEMALLOC(1,Paso_Solver_AMG); |
| 89 |
out->AMG_of_Schur=NULL; |
| 90 |
out->inv_A_FF=NULL; |
| 91 |
out->A_FF_pivot=NULL; |
| 92 |
out->A_FC=NULL; |
| 93 |
out->A_CF=NULL; |
| 94 |
out->rows_in_F=NULL; |
| 95 |
out->rows_in_C=NULL; |
| 96 |
out->mask_F=NULL; |
| 97 |
out->mask_C=NULL; |
| 98 |
out->x_F=NULL; |
| 99 |
out->b_F=NULL; |
| 100 |
out->x_C=NULL; |
| 101 |
out->b_C=NULL; |
| 102 |
out->A=Paso_SparseMatrix_getReference(A_p); |
| 103 |
out->level=level; |
| 104 |
|
| 105 |
if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) { |
| 106 |
/* identify independend set of rows/columns */ |
| 107 |
time0=Paso_timer(); |
| 108 |
#pragma omp parallel for private(i) schedule(static) |
| 109 |
for (i=0;i<n;++i) mis_marker[i]=-1; |
| 110 |
/*Paso_Pattern_RS(A_p,mis_marker,0.25);*/ |
| 111 |
Paso_Pattern_coup(A_p,mis_marker); |
| 112 |
time2=Paso_timer()-time0; |
| 113 |
if (Paso_noError()) { |
| 114 |
#pragma omp parallel for private(i) schedule(static) |
| 115 |
for (i = 0; i < n; ++i) counter[i]=mis_marker[i]; |
| 116 |
out->n=n; |
| 117 |
out->n_block=n_block; |
| 118 |
out->n_F=Paso_Util_cumsum(n,counter); |
| 119 |
out->mask_F=MEMALLOC(n,index_t); |
| 120 |
out->rows_in_F=MEMALLOC(out->n_F,index_t); |
| 121 |
out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double); |
| 122 |
out->A_FF_pivot=NULL; /* later use for block size>3 */ |
| 123 |
if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->inv_A_FF) || Paso_checkPtr(out->rows_in_F) ) ) { |
| 124 |
#pragma omp parallel |
| 125 |
{ |
| 126 |
/* creates an index for F from mask */ |
| 127 |
#pragma omp for private(i) schedule(static) |
| 128 |
for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1; |
| 129 |
#pragma omp for private(i) schedule(static) |
| 130 |
for (i = 0; i < n; ++i) { |
| 131 |
if (mis_marker[i]) { |
| 132 |
out->rows_in_F[counter[i]]=i; |
| 133 |
out->mask_F[i]=counter[i]; |
| 134 |
} else { |
| 135 |
out->mask_F[i]=-1; |
| 136 |
} |
| 137 |
} |
| 138 |
/* Compute row-sum for getting rs(A_FF)*/ |
| 139 |
#pragma omp for private(i,iPtr) schedule(static) |
| 140 |
for (i = 0; i < out->n_F; ++i) { |
| 141 |
out->inv_A_FF[i]=0; |
| 142 |
for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) { |
| 143 |
out->inv_A_FF[i]+=A_p->val[iPtr]; |
| 144 |
} |
| 145 |
} |
| 146 |
|
| 147 |
/* fprintf(stderr,"\n MATRIX START\n"); |
| 148 |
j=0; |
| 149 |
for (i = 0; i < n; ++i) { |
| 150 |
for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) { |
| 151 |
fprintf(stderr,"A[%d,%d]=%.2f ",i,A_p->pattern->index[iPtr],A_p->val[iPtr]); |
| 152 |
} |
| 153 |
if (mis_marker[i]) { |
| 154 |
fprintf(stderr,"ROW SUM %.2f\n ",out->inv_A_FF[j]); |
| 155 |
j++; |
| 156 |
} |
| 157 |
else { |
| 158 |
fprintf(stderr,"\n "); |
| 159 |
} |
| 160 |
} |
| 161 |
fprintf(stderr,"\n MATRIX END\n"); |
| 162 |
*/ |
| 163 |
#pragma omp for private(i, where_p,iPtr,index) schedule(static) |
| 164 |
for (i = 0; i < out->n_F; i++) { |
| 165 |
/* find main diagonal */ |
| 166 |
iPtr=A_p->pattern->ptr[out->rows_in_F[i]]; |
| 167 |
index=&(A_p->pattern->index[iPtr]); |
| 168 |
where_p=(index_t*)bsearch(&out->rows_in_F[i], |
| 169 |
index, |
| 170 |
A_p->pattern->ptr[out->rows_in_F[i] + 1]-A_p->pattern->ptr[out->rows_in_F[i]], |
| 171 |
sizeof(index_t), |
| 172 |
Paso_comparIndex); |
| 173 |
if (where_p==NULL) { |
| 174 |
Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing."); |
| 175 |
} else { |
| 176 |
iPtr+=(index_t)(where_p-index); |
| 177 |
/* get inverse of A_FF block: */ |
| 178 |
S=out->inv_A_FF[i]; |
| 179 |
if (ABS(A_p->val[iPtr])>0.) { |
| 180 |
if(ABS(S)>0.) |
| 181 |
out->inv_A_FF[i]=1./S; |
| 182 |
} else { |
| 183 |
Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getAMG: Break-down in AMG decomposition: non-regular main diagonal block."); |
| 184 |
} |
| 185 |
} |
| 186 |
} |
| 187 |
} /* end parallel region */ |
| 188 |
|
| 189 |
if( Paso_noError()) { |
| 190 |
/* if there are no nodes in the coarse level there is no more work to do */ |
| 191 |
out->n_C=n-out->n_F; |
| 192 |
/*if (out->n_C>11) {*/ |
| 193 |
/*fprintf(stderr,"SPARSITY in level %d \n", level);*/ |
| 194 |
if (level>0) { |
| 195 |
out->rows_in_C=MEMALLOC(out->n_C,index_t); |
| 196 |
out->mask_C=MEMALLOC(n,index_t); |
| 197 |
if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) { |
| 198 |
/* creates an index for C from mask */ |
| 199 |
#pragma omp parallel for private(i) schedule(static) |
| 200 |
for (i = 0; i < n; ++i) counter[i]=! mis_marker[i]; |
| 201 |
Paso_Util_cumsum(n,counter); |
| 202 |
#pragma omp parallel |
| 203 |
{ |
| 204 |
#pragma omp for private(i) schedule(static) |
| 205 |
for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1; |
| 206 |
#pragma omp for private(i) schedule(static) |
| 207 |
for (i = 0; i < n; ++i) { |
| 208 |
if (! mis_marker[i]) { |
| 209 |
out->rows_in_C[counter[i]]=i; |
| 210 |
out->mask_C[i]=counter[i]; |
| 211 |
} else { |
| 212 |
out->mask_C[i]=-1; |
| 213 |
} |
| 214 |
} |
| 215 |
} /* end parallel region */ |
| 216 |
/* get A_CF block: */ |
| 217 |
out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F); |
| 218 |
if (Paso_noError()) { |
| 219 |
/* get A_FC block: */ |
| 220 |
out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C); |
| 221 |
/* get A_CC block: */ |
| 222 |
if (Paso_noError()) { |
| 223 |
schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C); |
| 224 |
|
| 225 |
/*find the pattern of the schur complement with fill in*/ |
| 226 |
schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1); |
| 227 |
/* fprintf(stderr,"\nSHURE MATRIX START\n"); |
| 228 |
|
| 229 |
for (i = 0; i < schur->numRows; ++i) { |
| 230 |
for (iPtr=schur->pattern->ptr[i];iPtr<schur->pattern->ptr[i + 1]; ++iPtr) { |
| 231 |
fprintf(stderr,"A[%d,%d]=%.2f ",i,schur->pattern->index[iPtr],schur->val[iPtr]); |
| 232 |
} |
| 233 |
fprintf(stderr,"\n"); |
| 234 |
} |
| 235 |
fprintf(stderr,"MATRIX END\n"); |
| 236 |
|
| 237 |
fprintf(stderr,"\nA_CF MATRIX START\n"); |
| 238 |
for (i = 0; i < out->n_C; ++i) { |
| 239 |
for (iPtr=out->A_CF->pattern->ptr[i];iPtr<out->A_CF->pattern->ptr[i + 1]; ++iPtr) { |
| 240 |
fprintf(stderr,"A_CF[%d,%d]=%.2f ",i,out->A_CF->pattern->index[iPtr],out->A_CF->val[iPtr]); |
| 241 |
} |
| 242 |
fprintf(stderr,"\n"); |
| 243 |
} |
| 244 |
fprintf(stderr,"MATRIX END\n"); |
| 245 |
|
| 246 |
fprintf(stderr,"\nA_FC MATRIX START\n"); |
| 247 |
for (i = 0; i < out->n_F; ++i) { |
| 248 |
for (iPtr=out->A_FC->pattern->ptr[i];iPtr<out->A_FC->pattern->ptr[i + 1]; ++iPtr) { |
| 249 |
fprintf(stderr,"A_FC[%d,%d]=%.2f ",i,out->A_FC->pattern->index[iPtr],out->A_FC->val[iPtr]); |
| 250 |
} |
| 251 |
fprintf(stderr,"\n"); |
| 252 |
} |
| 253 |
fprintf(stderr,"MATRIX END\n"); |
| 254 |
|
| 255 |
mult=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern); |
| 256 |
fprintf(stderr,"\nPATTERN MATRIX START\n"); |
| 257 |
for (i = 0; i < mult->numOutput; ++i) { |
| 258 |
for (iPtr=mult->ptr[i];iPtr<mult->ptr[i + 1]; ++iPtr) { |
| 259 |
fprintf(stderr," P[%d,%d]=X ",i,mult->index[iPtr]); |
| 260 |
} |
| 261 |
fprintf(stderr,"\n"); |
| 262 |
} |
| 263 |
fprintf(stderr,"MATRIX END\n"); |
| 264 |
*/ |
| 265 |
/* copy values over*/ |
| 266 |
#pragma omp for private(i,iPtr,iPtr_s,j,j0) schedule(static) |
| 267 |
for (i = 0; i < schur_withFillIn->numRows; ++i) { |
| 268 |
for (iPtr=schur_withFillIn->pattern->ptr[i];iPtr<schur_withFillIn->pattern->ptr[i + 1]; ++iPtr) { |
| 269 |
j=schur_withFillIn->pattern->index[iPtr]; |
| 270 |
schur_withFillIn->val[iPtr]=0.; |
| 271 |
for (iPtr_s=schur->pattern->ptr[i];iPtr_s<schur->pattern->ptr[i + 1]; ++iPtr_s){ |
| 272 |
j0=schur->pattern->index[iPtr_s]; |
| 273 |
if (j==j0) { |
| 274 |
schur_withFillIn->val[iPtr]=schur->val[iPtr_s]; |
| 275 |
break; |
| 276 |
} |
| 277 |
} |
| 278 |
} |
| 279 |
} |
| 280 |
time0=Paso_timer()-time0; |
| 281 |
if (Paso_noError()) { |
| 282 |
time1=Paso_timer(); |
| 283 |
/* update A_CC block to get Schur complement and then apply AMG to it */ |
| 284 |
Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC); |
| 285 |
|
| 286 |
/* |
| 287 |
fprintf(stderr,"\n SHURE WITH FILL IN MATRIX START\n"); |
| 288 |
for (i = 0; i < schur_withFillIn->numRows; ++i) { |
| 289 |
for (iPtr=schur_withFillIn->pattern->ptr[i];iPtr<schur_withFillIn->pattern->ptr[i + 1]; ++iPtr) { |
| 290 |
fprintf(stderr,"A[%d,%d]=%.2f ",i,schur_withFillIn->pattern->index[iPtr],schur_withFillIn->val[iPtr]); |
| 291 |
} |
| 292 |
fprintf(stderr,"\n"); |
| 293 |
} |
| 294 |
fprintf(stderr,"MATRIX END\n"); |
| 295 |
*/ |
| 296 |
time1=Paso_timer()-time1; |
| 297 |
out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level-1); |
| 298 |
|
| 299 |
Paso_SparseMatrix_free(schur); |
| 300 |
/* Paso_SparseMatrix_free(schur_withFillIn);*/ |
| 301 |
} |
| 302 |
/* allocate work arrays for AMG application */ |
| 303 |
if (Paso_noError()) { |
| 304 |
out->x_F=MEMALLOC(n_block*out->n_F,double); |
| 305 |
out->b_F=MEMALLOC(n_block*out->n_F,double); |
| 306 |
out->x_C=MEMALLOC(n_block*out->n_C,double); |
| 307 |
out->b_C=MEMALLOC(n_block*out->n_C,double); |
| 308 |
|
| 309 |
if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) { |
| 310 |
#pragma omp parallel |
| 311 |
{ |
| 312 |
#pragma omp for private(i,k) schedule(static) |
| 313 |
for (i = 0; i < out->n_F; ++i) { |
| 314 |
for (k=0; k<n_block;++k) { |
| 315 |
out->x_F[i*n_block+k]=0.; |
| 316 |
out->b_F[i*n_block+k]=0.; |
| 317 |
} |
| 318 |
} |
| 319 |
#pragma omp for private(i,k) schedule(static) |
| 320 |
for (i = 0; i < out->n_C; ++i) { |
| 321 |
for (k=0; k<n_block;++k) { |
| 322 |
out->x_C[i*n_block+k]=0.; |
| 323 |
out->b_C[i*n_block+k]=0.; |
| 324 |
} |
| 325 |
} |
| 326 |
} /* end parallel region */ |
| 327 |
} |
| 328 |
} |
| 329 |
} |
| 330 |
} |
| 331 |
} |
| 332 |
} |
| 333 |
} |
| 334 |
} |
| 335 |
} |
| 336 |
} |
| 337 |
TMPMEMFREE(mis_marker); |
| 338 |
TMPMEMFREE(counter); |
| 339 |
if (Paso_noError()) { |
| 340 |
if (verbose) { |
| 341 |
printf("AMG: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F); |
| 342 |
if (out->n_C>0) { |
| 343 |
printf("timing: AMG: MIS/reordering/elemination : %e/%e/%e\n",time2,time0,time1); |
| 344 |
} else { |
| 345 |
printf("timing: AMG: MIS: %e\n",time2); |
| 346 |
} |
| 347 |
} |
| 348 |
return out; |
| 349 |
} else { |
| 350 |
Paso_Solver_AMG_free(out); |
| 351 |
return NULL; |
| 352 |
} |
| 353 |
} |
| 354 |
|
| 355 |
/**************************************************************/ |
| 356 |
|
| 357 |
/* apply AMG precondition b-> x |
| 358 |
|
| 359 |
in fact it solves |
| 360 |
|
| 361 |
[ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ] [ x_F ] = [b_F] |
| 362 |
[ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] [ x_C ] = [b_C] |
| 363 |
|
| 364 |
in the form |
| 365 |
|
| 366 |
b->[b_F,b_C] |
| 367 |
x_F=invA_FF*b_F |
| 368 |
b_C=b_C-A_CF*x_F |
| 369 |
x_C=AMG(b_C) |
| 370 |
b_F=b_F-A_FC*x_C |
| 371 |
x_F=invA_FF*b_F |
| 372 |
x<-[x_F,x_C] |
| 373 |
|
| 374 |
should be called within a parallel region |
| 375 |
barrier synconization should be performed to make sure that the input vector available |
| 376 |
|
| 377 |
*/ |
| 378 |
|
| 379 |
void Paso_Solver_solveAMG(Paso_Solver_AMG * amg, double * x, double * b) { |
| 380 |
dim_t i,k; |
| 381 |
dim_t n_block=amg->n_block; |
| 382 |
double *r=MEMALLOC(amg->n,double); |
| 383 |
Paso_Solver_GS* GS=NULL; |
| 384 |
|
| 385 |
if (amg->level==0) { |
| 386 |
/* x=invA_FF*b */ |
| 387 |
/*Paso_Solver_applyBlockDiagonalMatrix(n_block,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,x,b);*/ |
| 388 |
GS=Paso_Solver_getGS(amg->A,-1); |
| 389 |
Paso_Solver_solveGS(GS,x,b); |
| 390 |
Paso_Solver_GS_free(GS); |
| 391 |
} else { |
| 392 |
|
| 393 |
/* presmoothing */ |
| 394 |
GS=Paso_Solver_getGS(amg->A,-1); |
| 395 |
Paso_Solver_solveGS(GS,x,b); |
| 396 |
|
| 397 |
#pragma omp parallel for private(i) schedule(static) |
| 398 |
for (i=0;i<amg->n;++i) r[i]=b[i]; |
| 399 |
|
| 400 |
/*r=b-Ax*/ |
| 401 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r); |
| 402 |
/****************/ |
| 403 |
/* b->[b_F,b_C] */ |
| 404 |
if (n_block==1) { |
| 405 |
#pragma omp parallel for private(i) schedule(static) |
| 406 |
for (i=0;i<amg->n_F;++i) amg->b_F[i]=r[amg->rows_in_F[i]]; |
| 407 |
#pragma omp parallel for private(i) schedule(static) |
| 408 |
for (i=0;i<amg->n_C;++i) amg->b_C[i]=r[amg->rows_in_C[i]]; |
| 409 |
} else { |
| 410 |
#pragma omp parallel for private(i,k) schedule(static) |
| 411 |
for (i=0;i<amg->n_F;++i) |
| 412 |
for (k=0;k<n_block;k++) amg->b_F[amg->n_block*i+k]=r[n_block*amg->rows_in_F[i]+k]; |
| 413 |
#pragma omp parallel for private(i,k) schedule(static) |
| 414 |
for (i=0;i<amg->n_C;++i) |
| 415 |
for (k=0;k<n_block;k++) amg->b_C[amg->n_block*i+k]=r[n_block*amg->rows_in_C[i]+k]; |
| 416 |
} |
| 417 |
|
| 418 |
/* x_F=invA_FF*b_F */ |
| 419 |
Paso_Solver_applyBlockDiagonalMatrix(n_block,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F); |
| 420 |
|
| 421 |
/* b_C=b_C-A_CF*x_F */ |
| 422 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_CF,amg->x_F,1.,amg->b_C); |
| 423 |
|
| 424 |
/* x_C=AMG(b_C) */ |
| 425 |
Paso_Solver_solveAMG(amg->AMG_of_Schur,amg->x_C,amg->b_C); |
| 426 |
|
| 427 |
/* b_F=b_F-A_FC*x_C */ |
| 428 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_FC,amg->x_C,1.,amg->b_F); |
| 429 |
/* x_F=invA_FF*b_F */ |
| 430 |
Paso_Solver_applyBlockDiagonalMatrix(n_block,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F); |
| 431 |
/* x<-[x_F,x_C] */ |
| 432 |
|
| 433 |
if (n_block==1) { |
| 434 |
#pragma omp parallel for private(i) schedule(static) |
| 435 |
for (i=0;i<amg->n;++i) { |
| 436 |
if (amg->mask_C[i]>-1) { |
| 437 |
x[i]+=amg->x_C[amg->mask_C[i]]; |
| 438 |
} else { |
| 439 |
x[i]+=amg->x_F[amg->mask_F[i]]; |
| 440 |
} |
| 441 |
} |
| 442 |
} else { |
| 443 |
#pragma omp parallel for private(i,k) schedule(static) |
| 444 |
for (i=0;i<amg->n;++i) { |
| 445 |
if (amg->mask_C[i]>-1) { |
| 446 |
for (k=0;k<n_block;k++) x[n_block*i+k]+=amg->x_C[n_block*amg->mask_C[i]+k]; |
| 447 |
} else { |
| 448 |
for (k=0;k<n_block;k++) x[n_block*i+k]+=amg->x_F[n_block*amg->mask_F[i]+k]; |
| 449 |
} |
| 450 |
} |
| 451 |
} |
| 452 |
/* all done */ |
| 453 |
/*post smoothing*/ |
| 454 |
Paso_Solver_solveGS1(GS,x,b); |
| 455 |
|
| 456 |
Paso_Solver_GS_free(GS); |
| 457 |
} |
| 458 |
MEMFREE(r); |
| 459 |
return; |
| 460 |
} |
| 461 |
|
| 462 |
/* |
| 463 |
* $Log$ |
| 464 |
* |
| 465 |
*/ |
Parent Directory
| 
Revision Log