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

Revision 3323 - (show annotations)
Thu Oct 28 09:53:46 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 10593 byte(s)
```some openmp fizes.
```
 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 nessecary 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 ar 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 gives the new column number in P. S defines the strong connections. 39 40 If row i is in C, then P[i,j]=counter_C[i] if counter_C[i]>=0 and 0 otherwise 41 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 42 43 and a[i]= 44 alpha[i] = sum_s min(A[i,s],0)/(sum_{s in S and C} min(A[i,s],0)) A[i,k]<0 45 or if 46 beta[i] = sum_s max(A[i,s],0)/(sum_{s in S and C} max(A[i,s],0)) A[i,k]>0 47 48 49 */ 50 51 Paso_SparseMatrix* Paso_Preconditioner_AMG_getDirectProlongation(const Paso_SparseMatrix* A_p, 52 const dim_t* degree, const index_t* S, 53 const dim_t n_C, const index_t* counter_C) 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=A_p->pattern->ptr[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=A_p->pattern->ptr[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 (n_block == 1) { 120 Paso_Preconditioner_AMG_setDirectProlongation(out, A_p, counter_C); 121 } else { 122 Paso_Preconditioner_AMG_setDirectProlongation_Block(out, A_p, counter_C); 123 } 124 } 125 Paso_Pattern_free(outpattern); 126 if (Esys_noError()) { 127 return out; 128 } else { 129 Paso_SparseMatrix_free(out); 130 return NULL; 131 } 132 133 } 134 135 void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SparseMatrix* P_p, 136 const Paso_SparseMatrix* A_p, 137 const index_t *counter_C) { 138 dim_t i; 139 const dim_t n =A_p->numRows; 140 register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii; 141 register index_t iPtr, j, offset; 142 index_t *where_p, *start_p; 143 144 #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 ) schedule(static) 145 for (i=0;i=0) { 147 offset = P_p->pattern->ptr[i]; 148 P_p->val[offset]=1.; /* i is a C row */ 149 } else { 150 /* if i is an F row we first calculate alpha and beta: */ 151 sum_all_neg=0; /* sum of all negative values in row i of A */ 152 sum_all_pos=0; /* sum of all positive values in row i of A */ 153 sum_strong_neg=0; /* sum of all negative values A_ij where j is in C and strongly connected to i*/ 154 sum_strong_pos=0; /* sum of all positive values A_ij where j is in C and strongly connected to i*/ 155 A_ii=0; 156 for (iPtr=A_p->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 157 j=A_p->pattern->index[iPtr]; 158 A_ij=A_p->val[iPtr]; 159 if(j==i) { 160 A_ii=A_ij; 161 } else { 162 163 if(A_ij< 0) { 164 sum_all_neg+=A_ij; 165 } else { 166 sum_all_pos+=A_ij; 167 } 168 169 if (counter_C[j]>=0) { 170 /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */ 171 start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]); 172 where_p=(index_t*)bsearch(&(counter_C[j]), start_p, 173 P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i], 174 sizeof(index_t), 175 Paso_comparIndex); 176 if (! (where_p == NULL) ) { /* yes i stronly connect with j */ 177 offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p); 178 P_p->val[offset]=A_ij; /* will be modified later */ 179 if (A_ij< 0) { 180 sum_strong_neg+=A_ij; 181 } else { 182 sum_strong_pos+=A_ij; 183 } 184 } 185 } 186 187 } 188 } 189 if(sum_strong_neg<0) { 190 alpha= sum_all_neg/sum_strong_neg; 191 } else { 192 alpha=0; 193 } 194 if(sum_strong_pos>0) { 195 beta= sum_all_pos/sum_strong_pos; 196 } else { 197 beta=0; 198 A_ii+=sum_all_pos; 199 } 200 if ( A_ii > 0.) { 201 alpha*=(-1./A_ii); 202 beta*=(-1./A_ii); 203 } 204 205 for (iPtr=P_p->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 206 A_ij=P_p->val[iPtr]; 207 if (A_ij > 0 ) { 208 P_p->val[iPtr]*=beta; 209 } else { 210 P_p->val[iPtr]*=alpha; 211 } 212 } 213 } 214 } 215 } 216 217 void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SparseMatrix* P_p, 218 const Paso_SparseMatrix* A_p, 219 const index_t *counter_C) { 220 dim_t i; 221 const dim_t n =A_p->numRows; 222 const dim_t row_block=A_p->row_block_size; 223 const dim_t A_block = A_p->block_size; 224 double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii; 225 register double A_ij; 226 register index_t iPtr, j, offset, ib; 227 index_t *where_p, *start_p; 228 229 #pragma omp parallel private(ib, 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 ) 230 { 231 sum_all_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values in row i of A */ 232 sum_all_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values in row i of A */ 233 sum_strong_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values A_ij where j is in C and strongly connected to i*/ 234 sum_strong_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values A_ij where j is in C and strongly connected to i*/ 235 alpha=TMPMEMALLOC(row_block, double); 236 beta=TMPMEMALLOC(row_block, double); 237 A_ii=TMPMEMALLOC(row_block, double); 238 239 #pragma omp for schedule(static) 240 for (i=0;i=0) { 242 offset = P_p->pattern->ptr[i]; 243 for (ib =0; ibval[row_block*offset+ib]=1.; /* i is a C row */ 244 } else { 245 /* if i is an F row we first calculate alpha and beta: */ 246 for (ib =0; ibpattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 254 j=A_p->pattern->index[iPtr]; 255 if(j==i) { 256 for (ib =0; ibval[A_block*iPtr+ib+row_block*ib]; 257 } else { 258 for (ib =0; ibval[A_block*iPtr+ib+row_block*ib]; 260 if(A_ij< 0) { 261 sum_all_neg[ib]+=A_ij; 262 } else { 263 sum_all_pos[ib]+=A_ij; 264 } 265 } 266 267 if (counter_C[j]>=0) { 268 /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */ 269 start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]); 270 where_p=(index_t*)bsearch(&(counter_C[j]), start_p, 271 P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i], 272 sizeof(index_t), 273 Paso_comparIndex); 274 if (! (where_p == NULL) ) { /* yes i stronly connect with j */ 275 offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p); 276 for (ib =0; ibval[A_block*iPtr+ib+row_block*ib]; 278 P_p->val[row_block*offset+ib]=A_ij; /* will be modified later */ 279 if (A_ij< 0) { 280 sum_strong_neg[ib]+=A_ij; 281 } else { 282 sum_strong_pos[ib]+=A_ij; 283 } 284 } 285 } 286 } 287 288 } 289 } 290 for (ib =0; ib0) { 297 beta[ib]= sum_all_pos[ib]/sum_strong_pos[ib]; 298 } else { 299 beta[ib]=0; 300 A_ii[ib]+=sum_all_pos[ib]; 301 } 302 if ( A_ii[ib] > 0.) { 303 alpha[ib]*=(-1./A_ii[ib]); 304 beta[ib]*=(-1./A_ii[ib]); 305 } 306 } 307 308 for (iPtr=P_p->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 309 for (ib =0; ibval[row_block*iPtr+ib]; 311 if (A_ij > 0 ) { 312 P_p->val[row_block*iPtr+ib]*=beta[ib]; 313 } else { 314 P_p->val[row_block*iPtr+ib]*=alpha[ib]; 315 } 316 } 317 } 318 } 319 }/* end i loop */ 320 TMPMEMFREE(sum_all_neg); 321 TMPMEMFREE(sum_all_pos); 322 TMPMEMFREE(sum_strong_neg); 323 TMPMEMFREE(sum_strong_pos); 324 TMPMEMFREE(alpha); 325 TMPMEMFREE(beta); 326 TMPMEMFREE(A_ii); 327 } /* end parallel region */ 328 }