# Contents of /trunk/paso/src/AMG_Prolongation.c

Revision 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 23740 byte(s)
```Assorted spelling/comment fixes in paso.

```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2010 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: defines AMG prolongation */ 18 19 /**************************************************************/ 20 21 /* Author: Artak Amirbekyan, artak@uq.edu.au, l.gross@uq.edu.au */ 22 23 /**************************************************************/ 24 25 #include "Paso.h" 26 #include "SparseMatrix.h" 27 #include "PasoUtil.h" 28 #include "Preconditioner.h" 29 30 /************************************************************** 31 32 Methods necessary for AMG preconditioner 33 34 construct n x n_C the prolongation matrix P from A_p. 35 36 The columns in A_p to be considered are marked by counter_C[n] where 37 an unknown i to be considered in P is marked by 0<= counter_C[i] < n_C 38 and counter_C[i] gives the new column number in P. 39 S defines the strong connections. 40 41 The pattern of P is formed as follows: 42 43 If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise. 44 If row i is not C, then P[i,j] <> 0 if counter_C[k]==j (k in C) and (i,k) is a strong connection. 45 46 Two settings for P are implemented (see below) 47 48 */ 49 50 51 Paso_SparseMatrix* Paso_Preconditioner_AMG_getProlongation(Paso_SparseMatrix* A_p, 52 const index_t* offset_S, const dim_t* degree_S, const index_t* S, 53 const dim_t n_C, const index_t* counter_C, const index_t interpolation_method) 54 { 55 Paso_SparseMatrix* out=NULL; 56 Paso_Pattern *outpattern=NULL; 57 const dim_t n_block=A_p->row_block_size; 58 index_t *ptr=NULL, *index=NULL,j, iptr; 59 const dim_t n =A_p->numRows; 60 dim_t i,p,z, len_P; 61 62 ptr=MEMALLOC(n+1,index_t); 63 if (! Esys_checkPtr(ptr)) { 64 65 66 /* count the number of entries per row in the Prolongation matrix :*/ 67 68 #pragma omp parallel for private(i,z,iptr,j,p) schedule(static) 69 for (i=0;i=0) { 71 z=1; /* i is a C unknown */ 72 } else { 73 z=0; 74 iptr=offset_S[i]; 75 for (p=0; p=0) z++; /* and is in C */ 78 } 79 } 80 ptr[i]=z; 81 } 82 len_P=Paso_Util_cumsum(n,ptr); 83 ptr[n]=len_P; 84 85 /* allocate and create index vector for prolongation: */ 86 index=MEMALLOC(len_P,index_t); 87 88 if (! Esys_checkPtr(index)) { 89 #pragma omp parallel for private(i,z,iptr,j,p) schedule(static) 90 for (i=0;i=0) { 92 index[ptr[i]]=counter_C[i]; /* i is a C unknown */ 93 } else { 94 z=0; 95 iptr=offset_S[i]; 96 for (p=0; p=0) { /* and is in C */ 99 index[ptr[i]+z]=counter_C[j]; 100 z++; /* and is in C */ 101 } 102 } 103 } 104 } 105 } 106 } 107 if (Esys_noError()) { 108 outpattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,n,n_C,ptr,index); 109 } else { 110 MEMFREE(ptr); 111 MEMFREE(index); 112 } 113 /* now we need to create a matrix and fill it */ 114 if (Esys_noError()) { 115 out=Paso_SparseMatrix_alloc(MATRIX_FORMAT_DIAGONAL_BLOCK,outpattern,n_block,n_block,FALSE); 116 } 117 118 if (Esys_noError()) { 119 if ( (interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING) || ( interpolation_method == PASO_CLASSIC_INTERPOLATION) ) { 120 if (n_block == 1) { 121 Paso_Preconditioner_AMG_setClassicProlongation(out, A_p, offset_S, degree_S, S, counter_C); 122 } else { 123 Paso_Preconditioner_AMG_setClassicProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C); 124 } 125 } else { 126 if (n_block == 1) { 127 Paso_Preconditioner_AMG_setDirectProlongation(out, A_p, counter_C); 128 } else { 129 Paso_Preconditioner_AMG_setDirectProlongation_Block(out, A_p, counter_C); 130 } 131 } 132 } 133 Paso_Pattern_free(outpattern); 134 if (Esys_noError()) { 135 return out; 136 } else { 137 Paso_SparseMatrix_free(out); 138 return NULL; 139 } 140 141 } 142 143 /* 144 Direct Prolongation: 145 ------------------- 146 147 If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise. 148 If row i is not C, then P[i,j] = - a[i] * A[i,k]/A[i,i] with j=counter_C[k]>=0 and k in S 149 150 and a[i]= 151 alpha[i] = sum_s min(A[i,s],0)/(sum_{s in S and C} min(A[i,s],0)) A[i,k]<0 152 or if 153 beta[i] = sum_s max(A[i,s],0)/(sum_{s in S and C} max(A[i,s],0)) A[i,k]>0 154 155 156 */ 157 158 void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SparseMatrix* P_p, 159 const Paso_SparseMatrix* A_p, 160 const index_t *counter_C) { 161 dim_t i; 162 const dim_t n =A_p->numRows; 163 register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii, rtmp; 164 register index_t iPtr, j, offset; 165 index_t *where_p, *start_p; 166 167 #pragma omp parallel for private(A_ii, offset, where_p, start_p, i, alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos,iPtr,j, A_ij , rtmp) schedule(static) 168 for (i=0;i=0) { 170 offset = P_p->pattern->ptr[i]; 171 P_p->val[offset]=1.; /* i is a C row */ 172 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) { 173 /* if i is an F row we first calculate alpha and beta: */ 174 sum_all_neg=0; /* sum of all negative values in row i of A */ 175 sum_all_pos=0; /* sum of all positive values in row i of A */ 176 sum_strong_neg=0; /* sum of all negative values A_ij where j is in C and strongly connected to i*/ 177 sum_strong_pos=0; /* sum of all positive values A_ij where j is in C and strongly connected to i*/ 178 A_ii=0; 179 for (iPtr=A_p->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 180 j=A_p->pattern->index[iPtr]; 181 A_ij=A_p->val[iPtr]; 182 if(j==i) { 183 A_ii=A_ij; 184 } else { 185 186 if(A_ij< 0) { 187 sum_all_neg+=A_ij; 188 } else { 189 sum_all_pos+=A_ij; 190 } 191 192 if (counter_C[j]>=0) { 193 /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */ 194 start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]); 195 where_p=(index_t*)bsearch(&(counter_C[j]), start_p, 196 P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i], 197 sizeof(index_t), 198 Paso_comparIndex); 199 if (! (where_p == NULL) ) { /* yes i strongly connected with j */ 200 offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p); 201 P_p->val[offset]=A_ij; /* will be modified later */ 202 if (A_ij< 0) { 203 sum_strong_neg+=A_ij; 204 } else { 205 sum_strong_pos+=A_ij; 206 } 207 } 208 } 209 210 } 211 } 212 if(sum_strong_neg<0) { 213 alpha= sum_all_neg/sum_strong_neg; 214 } else { 215 alpha=0; 216 } 217 if(sum_strong_pos>0) { 218 beta= sum_all_pos/sum_strong_pos; 219 } else { 220 beta=0; 221 A_ii+=sum_all_pos; 222 } 223 if ( A_ii > 0.) { 224 rtmp=(-1.)/A_ii; 225 alpha*=rtmp; 226 beta*=rtmp; 227 } 228 for (iPtr=P_p->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 229 A_ij=P_p->val[iPtr]; 230 if (A_ij > 0 ) { 231 P_p->val[iPtr]*=beta; 232 } else { 233 P_p->val[iPtr]*=alpha; 234 } 235 } 236 } 237 } 238 } 239 240 void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SparseMatrix* P_p, 241 const Paso_SparseMatrix* A_p, 242 const index_t *counter_C) { 243 dim_t i; 244 const dim_t n =A_p->numRows; 245 const dim_t row_block=A_p->row_block_size; 246 const dim_t A_block = A_p->block_size; 247 double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii; 248 register double A_ij, rtmp; 249 register index_t iPtr, j, offset, ib; 250 index_t *where_p, *start_p; 251 252 #pragma omp parallel private(ib, rtmp, A_ii, offset, where_p, start_p, i, alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos,iPtr,j, A_ij ) 253 { 254 sum_all_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values in row i of A */ 255 sum_all_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values in row i of A */ 256 sum_strong_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values A_ij where j is in C and strongly connected to i*/ 257 sum_strong_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values A_ij where j is in C and strongly connected to i*/ 258 alpha=TMPMEMALLOC(row_block, double); 259 beta=TMPMEMALLOC(row_block, double); 260 A_ii=TMPMEMALLOC(row_block, double); 261 262 #pragma omp for schedule(static) 263 for (i=0;i=0) { 265 offset = P_p->pattern->ptr[i]; 266 for (ib =0; ibval[row_block*offset+ib]=1.; /* i is a C row */ 267 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) { 268 /* if i is an F row we first calculate alpha and beta: */ 269 for (ib =0; ibpattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 277 j=A_p->pattern->index[iPtr]; 278 if(j==i) { 279 for (ib =0; ibval[A_block*iPtr+ib+row_block*ib]; 280 } else { 281 for (ib =0; ibval[A_block*iPtr+ib+row_block*ib]; 283 if(A_ij< 0) { 284 sum_all_neg[ib]+=A_ij; 285 } else { 286 sum_all_pos[ib]+=A_ij; 287 } 288 } 289 290 if (counter_C[j]>=0) { 291 /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */ 292 start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]); 293 where_p=(index_t*)bsearch(&(counter_C[j]), start_p, 294 P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i], 295 sizeof(index_t), 296 Paso_comparIndex); 297 if (! (where_p == NULL) ) { /* yes i strongly connected with j */ 298 offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p); 299 for (ib =0; ibval[A_block*iPtr+ib+row_block*ib]; 301 P_p->val[row_block*offset+ib]=A_ij; /* will be modified later */ 302 if (A_ij< 0) { 303 sum_strong_neg[ib]+=A_ij; 304 } else { 305 sum_strong_pos[ib]+=A_ij; 306 } 307 } 308 } 309 } 310 311 } 312 } 313 for (ib =0; ib0) { 320 beta[ib]= sum_all_pos[ib]/sum_strong_pos[ib]; 321 } else { 322 beta[ib]=0; 323 A_ii[ib]+=sum_all_pos[ib]; 324 } 325 if ( A_ii[ib] > 0.) { 326 rtmp=(-1./A_ii[ib]); 327 alpha[ib]*=rtmp; 328 beta[ib]*=rtmp; 329 } 330 } 331 332 for (iPtr=P_p->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 333 for (ib =0; ibval[row_block*iPtr+ib]; 335 if (A_ij > 0 ) { 336 P_p->val[row_block*iPtr+ib]*=beta[ib]; 337 } else { 338 P_p->val[row_block*iPtr+ib]*=alpha[ib]; 339 } 340 } 341 } 342 } 343 }/* end i loop */ 344 TMPMEMFREE(sum_all_neg); 345 TMPMEMFREE(sum_all_pos); 346 TMPMEMFREE(sum_strong_neg); 347 TMPMEMFREE(sum_strong_pos); 348 TMPMEMFREE(alpha); 349 TMPMEMFREE(beta); 350 TMPMEMFREE(A_ii); 351 } /* end parallel region */ 352 } 353 354 /* 355 Classic Prolongation: 356 ------------------- 357 358 If row i is in C (counter_C[i]>=0), then P[i,j]=1 if j==counter_C[i] or 0 otherwise. 359 If row i is not C, then P[i,j] = - 1/a[i] * ( A[i,k] + sum_{l} A[i,l]*A+[l,k]/B[i,k]) 360 where the summation over l is considering columns which are strongly connected 361 to i (l in S[i]) and not in C (counter_C[l]<0) and 362 363 B[i,k]=sum_{m in S_i and in C} A+[k,m] 364 a[i]=A[i,i]+sum{l not strongly connected to i} A[i,l] 365 366 A+[i,k]=A[i,k] if sign(A[i,k])==sign(A[i,i]) or 0 otherwise 367 368 369 */ 370 void Paso_Preconditioner_AMG_setClassicProlongation(Paso_SparseMatrix* P_p, 371 Paso_SparseMatrix* A_p, 372 const index_t* offset_S, const dim_t* degree_S, const index_t* S, 373 const index_t *counter_C) { 374 dim_t i, q; 375 const dim_t n =A_p->numRows; 376 double *D_s=NULL; 377 index_t *D_s_offset=NULL, iPtr, iPtr_j; 378 const dim_t ll = Paso_Util_iMax(n, degree_S); 379 const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p); 380 381 382 #pragma omp parallel private(D_s, D_s_offset, iPtr, q, iPtr_j) 383 { 384 D_s=TMPMEMALLOC(ll,double); 385 D_s_offset=TMPMEMALLOC(ll,index_t); 386 387 388 #pragma omp for private(i) schedule(static) 389 for (i=0;i=0) { 391 P_p->val[P_p->pattern->ptr[i]]=1.; /* i is a C row */ 392 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) { 393 const index_t *start_s = &(S[offset_S[i]]); 394 const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]); 395 const dim_t degree_P_i = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i]; 396 /* this loop sums up the weak connections in A and creates a list of the strong connected columns 397 which are not in C (=no interpolation nodes) */ 398 const double A_ii = A_p->val[ptr_main_A[i]]; 399 double a=A_ii; 400 401 for (iPtr=A_p->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 402 const index_t j=A_p->pattern->index[iPtr]; 403 const double A_ij=A_p->val[iPtr]; 404 if ( (i!=j) && (degree_S[j]>0) ) { 405 /* is (i,j) a strong connection ?*/ 406 const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex); 407 if (where_s == NULL) { /* weak connections are accumulated */ 408 a+=A_ij; 409 } else { /* yes i strongly connected with j */ 410 if (counter_C[j]>=0) { /* j is an interpolation point : add A_ij into P */ 411 const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex); 412 if (where_p == NULL) { 413 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation: interpolation point is missing."); 414 } else { 415 const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p); 416 P_p->val[offset]+=A_ij; 417 } 418 } else { /* j is not an interpolation point */ 419 /* find all interpolation points m of k */ 420 const index_t *start_p_j = &(P_p->pattern->index[P_p->pattern->ptr[i]]); 421 const dim_t degree_P_j = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i]; 422 double s=0.; 423 dim_t len_D_s=0; 424 for (iPtr_j=A_p->pattern->ptr[j];iPtr_jpattern->ptr[j + 1]; ++iPtr_j) { 425 const double A_jm=A_p->val[iPtr_j]; 426 const index_t m=A_p->pattern->index[iPtr_j]; 427 /* is m an interpolation point ? */ 428 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_j,degree_P_j, sizeof(index_t), Paso_comparIndex); 429 if (! (where_p_m==NULL)) { 430 const index_t offset_m = P_p->pattern->ptr[i]+ (index_t)(where_p_m-start_p_j); 431 if (! SAMESIGN(A_ii,A_jm)) { 432 D_s[len_D_s]=A_jm; 433 } else { 434 D_s[len_D_s]=0.; 435 } 436 D_s_offset[len_D_s]=offset_m; 437 len_D_s++; 438 } 439 } 440 for (q=0;q0) { 442 s=A_ij/s; 443 for (q=0;qval[D_s_offset[q]]+=s*D_s[q]; 445 } 446 } else { 447 a+=A_ij; 448 } 449 } 450 } 451 } 452 } /* i has been processed, now we need to do some rescaling */ 453 if (ABS(a)>0.) { 454 a=-1./a; 455 for (iPtr=P_p->pattern->ptr[i]; iPtrpattern->ptr[i + 1]; ++iPtr) { 456 P_p->val[iPtr]*=a; 457 } 458 } 459 } 460 } /* end of row i loop */ 461 TMPMEMFREE(D_s); 462 TMPMEMFREE(D_s_offset); 463 } /* end of parallel region */ 464 } 465 466 void Paso_Preconditioner_AMG_setClassicProlongation_Block(Paso_SparseMatrix* P_p, 467 Paso_SparseMatrix* A_p, 468 const index_t* offset_S, const dim_t* degree_S, const index_t* S, 469 const index_t *counter_C) { 470 dim_t i, q, ib; 471 const dim_t row_block=A_p->row_block_size; 472 const dim_t A_block = A_p->block_size; 473 const dim_t n =A_p->numRows; 474 double *D_s=NULL; 475 index_t *D_s_offset=NULL, iPtr, iPtr_j; 476 const dim_t ll = Paso_Util_iMax(n, degree_S); 477 const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p); 478 479 480 #pragma omp parallel private(D_s, D_s_offset, iPtr, q, iPtr_j,ib) 481 { 482 double *a=TMPMEMALLOC(row_block, double); 483 D_s=TMPMEMALLOC(row_block*ll,double); 484 D_s_offset=TMPMEMALLOC(row_block*ll,index_t); 485 486 #pragma omp for private(i) schedule(static) 487 for (i=0;i=0) { 489 const index_t offset = P_p->pattern->ptr[i]; 490 for (ib =0; ibval[row_block*offset+ib]=1.; /* i is a C row */ 491 } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) { 492 const index_t *start_s = &(S[offset_S[i]]); 493 const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]); 494 const dim_t degree_P_i = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i]; 495 /* this loop sums up the weak connections in A and creates a list of the strong connected columns 496 which are not in C (=no interpolation nodes) */ 497 const double *A_ii = &(A_p->val[ptr_main_A[i]*A_block]); 498 for (ib=0; ibpattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 502 const index_t j=A_p->pattern->index[iPtr]; 503 const double* A_ij=&(A_p->val[iPtr*A_block]); 504 505 if ( (i!=j) && (degree_S[j]>0) ) { 506 /* is (i,j) a strong connection ?*/ 507 const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex); 508 if (where_s == NULL) { /* weak connections are accumulated */ 509 for (ib=0; ib=0) { /* j is an interpolation point : add A_ij into P */ 512 const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex); 513 if (where_p == NULL) { 514 Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation: interpolation point is missing."); 515 } else { 516 const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p); 517 for (ib=0; ibval[offset*row_block+ib] +=A_ij[(row_block+1)*ib]; 518 } 519 } else { /* j is not an interpolation point */ 520 /* find all interpolation points m of k */ 521 const index_t *start_p_j = &(P_p->pattern->index[P_p->pattern->ptr[i]]); 522 const dim_t degree_P_j = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i]; 523 dim_t len_D_s=0; 524 for (iPtr_j=A_p->pattern->ptr[j];iPtr_jpattern->ptr[j + 1]; ++iPtr_j) { 525 const double* A_jm=&(A_p->val[iPtr_j*A_block]); 526 const index_t m=A_p->pattern->index[iPtr_j]; 527 /* is m an interpolation point ? */ 528 const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_j,degree_P_j, sizeof(index_t), Paso_comparIndex); 529 if (! (where_p_m==NULL)) { 530 const index_t offset_m = P_p->pattern->ptr[i]+ (index_t)(where_p_m-start_p_j); 531 for (ib=0; ib0) { 547 s=A_ij[(row_block+1)*ib]/s; 548 for (q=0;qval[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib]; 550 } 551 } else { 552 a[ib]+=A_ij[(row_block+1)*ib]; 553 } 554 } 555 } 556 } 557 } 558 } /* i has been processed, now we need to do some rescaling */ 559 for (ib=0; ib0.) { 562 a2=-1./a2; 563 for (iPtr=P_p->pattern->ptr[i]; iPtrpattern->ptr[i + 1]; ++iPtr) { 564 P_p->val[iPtr*row_block+ib]*=a2; 565 } 566 } 567 } 568 } 569 } /* end of row i loop */ 570 TMPMEMFREE(D_s); 571 TMPMEMFREE(D_s_offset); 572 TMPMEMFREE(a); 573 } /* end of parallel region */ 574 }