Contents of /trunk/paso/src/SparseMatrix_AMGcomponents.c

Revision 2828 - (show annotations)
Tue Dec 22 01:24:40 2009 UTC (11 years, 5 months ago) by artak
File MIME type: text/plain
File size: 11099 byte(s)
```Smoother for AMG now can be selected from python. Now only Jacobi and Gauss-Seidel are available as smoothers.
```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2009 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: SparseMatrix_AMGcomponents */ 18 19 /**************************************************************/ 20 21 /* Author: Artak Amirbekyan, artak@uq.edu.au */ 22 23 /**************************************************************/ 24 25 #include "Paso.h" 26 #include "SparseMatrix.h" 27 #include "PasoUtil.h" 28 29 /************************************************************** 30 31 Methods nessecary for AMG preconditioner 32 33 */ 34 35 /* Prolongation matrix is constracted by W_FC as P<- [W_FC] 36 [I_CC] 37 */ 38 39 Paso_SparseMatrix* Paso_SparseMatrix_getProlongation(Paso_SparseMatrix* W, index_t* mis_marker){ 40 41 Paso_Pattern *outpattern=NULL; 42 Paso_SparseMatrix *out=NULL; 43 index_t iptr,wptr; 44 45 Paso_IndexList* index_list=NULL; 46 dim_t n=W->numRows+W->numCols; 47 dim_t i,j,k=0; 48 49 index_list=TMPMEMALLOC(n,Paso_IndexList); 50 if (! Paso_checkPtr(index_list)) { 51 #pragma omp parallel for private(i) schedule(static) 52 for(i=0;ipattern->ptr[k];iptrpattern->ptr[k+1]; ++iptr) { 62 j=W->pattern->index[iptr]; 63 Paso_IndexList_insertIndex(&(index_list[i]),j); 64 } 65 k++; 66 } 67 else { 68 Paso_IndexList_insertIndex(&(index_list[i]),i-k); 69 } 70 } 71 72 outpattern=Paso_IndexList_createPattern(0, n,index_list,0,W->numCols,0); 73 out=Paso_SparseMatrix_alloc(W->type,outpattern,1,1,TRUE); 74 75 k=0; 76 77 for (i=0;ipattern->ptr[k]; 80 for (iptr=out->pattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 81 out->val[iptr]=W->val[wptr]; 82 wptr++; 83 } 84 k++; 85 } 86 else { 87 iptr=out->pattern->ptr[i]; 88 out->val[iptr]=1; 89 } 90 } 91 92 /* clean up */ 93 if (index_list!=NULL) { 94 #pragma omp parallel for private(i) schedule(static) 95 for(i=0;inumCols; 112 dim_t F=P->numRows-C; 113 dim_t n=C+F; 114 dim_t i,j,k=0; 115 index_t iptr,jptr; 116 117 index_list=TMPMEMALLOC(C,Paso_IndexList); 118 if (! Paso_checkPtr(index_list)) { 119 #pragma omp parallel for private(i) schedule(static) 120 for(i=0;ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 129 j=P->pattern->index[iptr]; 130 Paso_IndexList_insertIndex(&(index_list[j]),i); 131 } 132 } 133 134 outpattern=Paso_IndexList_createPattern(0, C,index_list,0,C+F,0); 135 out=Paso_SparseMatrix_alloc(P->type,outpattern,1,1,TRUE); 136 137 for (i=0;inumRows;++i) { 138 for (iptr=out->pattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 139 j=out->pattern->index[iptr]; 140 /*This for later will be replaced by bsearch!!*/ 141 for (jptr=P->pattern->ptr[j];jptrpattern->ptr[j+1]; ++jptr) { 142 k=P->pattern->index[jptr]; 143 if(k==i) { 144 out->val[iptr]=P->val[jptr]; 145 } 146 } 147 } 148 } 149 150 /* clean up */ 151 if (index_list!=NULL) { 152 #pragma omp parallel for private(i) schedule(static) 153 for(i=0;inumRows; 169 dim_t F=W_FC->numRows; 170 dim_t i,j,k; 171 172 index_t iPtr,dptr=0; 173 /*index_t *index, *where_p;*/ 174 175 alpha=TMPMEMALLOC(F,double); 176 beta=TMPMEMALLOC(F,double); 177 178 179 k=0; 180 for (i = 0; i < n; ++i) { 181 if(mis_marker[i]) { 182 alpha[k]=0; 183 beta[k]=0; 184 sum_all_neg=0; 185 sum_all_pos=0; 186 sum_strong_neg=0; 187 sum_strong_pos=0; 188 /*compute beta_i=sum_{j\inN_i}a^{+}_ij and alpha_i=sum_{j\inN_i}a^{-}_ij and denominators for both alpha_i and beta_i*/ 189 for (iPtr=A->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 190 j=A->pattern->index[iPtr]; 191 if(j!=i) { 192 if(A->val[iPtr]<0) { 193 sum_all_neg+=A->val[iPtr]; 194 if(!mis_marker[j]) { 195 sum_strong_neg+=A->val[iPtr]; 196 } 197 } 198 else if(A->val[iPtr]>0) { 199 sum_all_pos+=A->val[iPtr]; 200 if(!mis_marker[j]) { 201 sum_strong_pos+=A->val[iPtr]; 202 } 203 } 204 } 205 else { 206 a_ii=A->val[iPtr]; 207 dptr=iPtr; 208 } 209 } 210 211 if(sum_strong_neg!=0) { 212 alpha[k]=sum_all_neg/(sum_strong_neg); 213 } 214 if(sum_strong_pos!=0) { 215 beta[k]=sum_all_pos/(sum_strong_pos); 216 } 217 else { 218 a_ii+=sum_all_pos; 219 beta[k]=0; 220 } 221 alpha[k]=alpha[k]/(a_ii); 222 beta[k]=beta[k]/(a_ii); 223 k++; 224 /*printf("Got in row=%d, alpha[%d]=%e, beta[%d]=%e, a_den=%e, b_den=%e \n",i,k-1,alpha[k-1],k-1,beta[k-1],alpha_den[k-1],beta_den[k-1]);*/ 225 } 226 } 227 #pragma omp parallel for private(i,iPtr) schedule(static) 228 for (i = 0; i < W_FC->numRows; ++i) { 229 for (iPtr=W_FC->pattern->ptr[i];iPtrpattern->ptr[i + 1]; ++iPtr) { 230 if(W_FC->val[iPtr]<0) { 231 W_FC->val[iPtr]=-alpha[i]*W_FC->val[iPtr]; 232 } 233 else if(W_FC->val[iPtr]>0) { 234 W_FC->val[iPtr]=-beta[i]*W_FC->val[iPtr]; 235 } 236 } 237 } 238 239 TMPMEMFREE(alpha); 240 TMPMEMFREE(beta); 241 } 242 243 244 /*A_c=R*A*P taking into account coarseneing */ 245 /*Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix *R, Paso_SparseMatrix *P) { 246 247 Paso_SparseMatrix *temp=NULL; 248 Paso_SparseMatrix *out=NULL; 249 250 temp=Paso_SparseMatrix_MatrixMatrix(A,P); 251 out=Paso_SparseMatrix_MatrixMatrix(R,temp); 252 253 Paso_SparseMatrix_free(temp); 254 255 return out; 256 } 257 */ 258 259 Paso_SparseMatrix* Paso_SparseMatrix_MatrixMatrix(Paso_SparseMatrix* A, Paso_SparseMatrix* B) { 260 index_t iptrA,iptrB,iptrC; 261 dim_t i,j=0,k,l; 262 Paso_Pattern* outpattern=NULL; 263 Paso_SparseMatrix *out=NULL; 264 double sum,b_lj=0; 265 266 double time0=0; 267 bool_t verbose=0; 268 269 time0=Paso_timer(); 270 271 outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,B->pattern); 272 out=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE); 273 274 time0=Paso_timer()-time0; 275 if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Pattern creation: %e\n",time0); 276 277 time0=Paso_timer(); 278 279 #pragma omp parallel for private(i,iptrC,j,sum,iptrA,k,b_lj,iptrB,l) schedule(static) 280 for(i = 0; i < out->numRows; i++) { 281 for(iptrC = out->pattern->ptr[i]; iptrC < out->pattern->ptr[i+1]; ++iptrC) { 282 j = out->pattern->index[iptrC]; 283 sum=0; 284 for(iptrA = A->pattern->ptr[i]; iptrA < A->pattern->ptr[i+1]; ++iptrA) { 285 k=A->pattern->index[iptrA]; 286 /*can be replaced by bsearch */ 287 b_lj=0; 288 for(iptrB = B->pattern->ptr[k]; iptrB < B->pattern->ptr[k+1]; ++iptrB) { 289 l=B->pattern->index[iptrB]; 290 if(l==j) { 291 b_lj=B->val[iptrB]; 292 break; 293 } 294 } 295 sum+=A->val[iptrA]*b_lj; 296 } 297 out->val[iptrC]=sum; 298 } 299 } 300 301 time0=Paso_timer()-time0; 302 if (verbose) fprintf(stdout,"timing: Paso_SparseMatrix_MatrixMatrix: Matrix multiplication: %e\n",time0); 303 304 Paso_Pattern_free(outpattern); 305 306 return out; 307 } 308 309 310 Paso_SparseMatrix* Paso_Solver_getCoarseMatrix(Paso_SparseMatrix* A,Paso_SparseMatrix *R,Paso_SparseMatrix *P) { 311 312 index_t iptrA_c,iptrR,iptrP,iptrA; 313 dim_t i,j,k,l,m; 314 double first_sum=0,second_sum=0,p_lj=0,a_kl=0,r_ik=0; 315 Paso_Pattern* temp=NULL; 316 Paso_Pattern* outpattern=NULL; 317 Paso_SparseMatrix *A_c=NULL; 318 319 double time0=0; 320 bool_t verbose=0; 321 322 time0=Paso_timer(); 323 324 temp=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,A->pattern,P->pattern); 325 outpattern=Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,R->pattern,temp); 326 A_c=Paso_SparseMatrix_alloc(A->type,outpattern,1,1, TRUE); 327 328 time0=Paso_timer()-time0; 329 if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Pattern creation: %e\n",time0); 330 331 /*a^c_ij=sum_k^n(r_ik)sum_l^n(a_kl*P_lj)*/ 332 333 time0=Paso_timer(); 334 335 #pragma omp parallel for private(i,iptrA_c,j,second_sum,iptrR,k,first_sum,p_lj,iptrP,m,a_kl,r_ik) schedule(static) 336 for(i = 0; i < A_c->numRows; i++) { 337 for(iptrA_c = A_c->pattern->ptr[i]; iptrA_c < A_c->pattern->ptr[i+1]; ++iptrA_c) { 338 j = A_c->pattern->index[iptrA_c]; 339 second_sum=0; 340 for(iptrR = R->pattern->ptr[i]; iptrR < R->pattern->ptr[i+1]; ++iptrR) { 341 k = R->pattern->index[iptrR]; 342 first_sum=0; 343 for(iptrA = A->pattern->ptr[k]; iptrA < A->pattern->ptr[k+1]; ++iptrA) { 344 l= A->pattern->index[iptrA]; 345 p_lj=0; 346 /*This loop can be replaced by bsearch */ 347 for(iptrP = P->pattern->ptr[l]; iptrP < P->pattern->ptr[l+1]; ++iptrP) { 348 m=P->pattern->index[iptrP]; 349 if(m==j) { 350 p_lj=P->val[iptrP]; 351 break; 352 } 353 } 354 a_kl=A->val[iptrA]; 355 first_sum+=a_kl*p_lj; 356 } 357 r_ik=R->val[iptrR]; 358 second_sum+=r_ik*first_sum; 359 } 360 A_c->val[iptrA_c]=second_sum; 361 } 362 } 363 time0=Paso_timer()-time0; 364 if (verbose) fprintf(stdout,"timing: Paso_Solver_getCoarseMatrix: Matrix multiplication: %e\n",time0); 365 366 Paso_Pattern_free(outpattern); 367 Paso_Pattern_free(temp); 368 369 return A_c; 370 } 371

 ViewVC Help Powered by ViewVC 1.1.26