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

Revision 1028 - (show annotations)
Wed Mar 14 00:15:24 2007 UTC (12 years, 11 months ago) by gross
File MIME type: text/plain
File size: 24375 byte(s)
```modifications to be compliant with _WIN32. The substitutes for asinh, acosh, atanh are still missing (erf will through an exception)
```
 1 /* \$Id\$ */ 2 3 4 /* 5 ******************************************************************************** 6 * Copyright 2006 by ACcESS MNRF * 7 * * 8 * http://www.access.edu.au * 9 * Primary Business: Queensland, Australia * 10 * Licensed under the Open Software License version 3.0 * 11 * http://www.opensource.org/licenses/osl-3.0.php * 12 ******************************************************************************** 13 */ 14 15 /* 16 * Purpose 17 * ======= 18 * 19 * GMRES solves the linear system A*x=b using the 20 * truncated and restered GMRES method with preconditioning. 21 * 22 * Convergence test: norm( b - A*x )< TOL. 23 * 24 * 25 * 26 * Arguments 27 * ========= 28 * 29 * r (input/output) double array, dimension n. 30 * On entry, residual of inital guess X 31 * 32 * x (input/output) double array, dimension n. 33 * On input, the initial guess. 34 * 35 * iter (input/output) int 36 * On input, the maximum num_iterations to be performed. 37 * On output, actual number of num_iterations performed. 38 * 39 * tolerance (input/output) DOUBLE PRECISION 40 * On input, the allowable convergence measure for 41 * norm( b - A*x ) 42 * On output, the final value of this measure. 43 * 44 * Length_of_recursion (input) gives the number of residual to be kept in orthogonalization process 45 * 46 * restart (input) If restart>0, iteration is resterted a after restart steps. 47 * 48 * INFO (output) int 49 * 50 * =SOLVER_NO_ERROR: Successful exit. num_iterated approximate solution returned. 51 * =SOLVER_MAXNUM_ITER_REACHED 52 * =SOLVER_INPUT_ERROR Illegal parameter: 53 * =SOLVER_BREAKDOWN: bad luck! 54 * =SOLVER_MEMORY_ERROR : no memory available 55 * 56 * ============================================================== 57 */ 58 59 #include "Common.h" 60 #include "SystemMatrix.h" 61 #include "Solver.h" 62 #ifdef _OPENMP 63 #include 64 #endif 65 66 err_t Paso_Solver_GMRES( 67 Paso_SystemMatrix * A, 68 double * r, 69 double * x, 70 dim_t *iter, 71 double * tolerance,dim_t Length_of_recursion,dim_t restart, 72 Paso_Performance* pp) { 73 74 /* Local variables */ 75 76 double *AP,**X_PRES,**R_PRES,**P_PRES; 77 double *P_PRES_dot_AP,*R_PRES_dot_P_PRES,*BREAKF,*ALPHA; 78 double P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3,P_PRES_dot_AP4,P_PRES_dot_AP5,P_PRES_dot_AP6,R_PRES_dot_P,abs_RP,breakf0; 79 double tol,Factor,sum_BREAKF,gamma,SC1,SC2,norm_of_residual,diff,L2_R,Norm_of_residual_global; 80 double *save_XPRES, *save_P_PRES, *save_R_PRES,save_R_PRES_dot_P_PRES; 81 dim_t maxit,Num_iter_global,num_iter_restart,num_iter; 82 dim_t i,z,order,n, Length_of_mem; 83 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE; 84 err_t Status=SOLVER_NO_ERROR; 85 86 87 /* adapt original routine parameters */ 88 89 n=A->num_cols * A-> col_block_size; 90 Length_of_mem=MAX(Length_of_recursion,0)+1; 91 92 93 /* Test the input parameters. */ 94 if (restart>0) restart=MAX(Length_of_recursion,restart); 95 if (n < 0 || Length_of_recursion<=0) { 96 return SOLVER_INPUT_ERROR; 97 } 98 99 /* allocate memory: */ 100 X_PRES=TMPMEMALLOC(Length_of_mem,double*); 101 R_PRES=TMPMEMALLOC(Length_of_mem,double*); 102 P_PRES=TMPMEMALLOC(Length_of_mem,double*); 103 P_PRES_dot_AP=TMPMEMALLOC(Length_of_mem,double); 104 R_PRES_dot_P_PRES=TMPMEMALLOC(Length_of_mem,double); 105 BREAKF=TMPMEMALLOC(Length_of_mem,double); 106 ALPHA=TMPMEMALLOC(Length_of_mem,double); 107 AP=TMPMEMALLOC(n,double); 108 if (AP==NULL || X_PRES ==NULL || R_PRES == NULL || P_PRES == NULL || 109 P_PRES_dot_AP == NULL || R_PRES_dot_P_PRES ==NULL || BREAKF == NULL || ALPHA == NULL) { 110 Status=SOLVER_MEMORY_ERROR; 111 } else { 112 for (i=0;i6) { 332 #pragma omp master 333 { 334 R_PRES_dot_P=ZERO; 335 P_PRES_dot_AP0=ZERO; 336 P_PRES_dot_AP1=ZERO; 337 P_PRES_dot_AP2=ZERO; 338 P_PRES_dot_AP3=ZERO; 339 P_PRES_dot_AP4=ZERO; 340 P_PRES_dot_AP5=ZERO; 341 P_PRES_dot_AP6=ZERO; 342 } 343 #pragma omp barrier 344 #pragma omp for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3,P_PRES_dot_AP4,P_PRES_dot_AP5,P_PRES_dot_AP6) schedule(static) 345 for (z=0;z ZERO) && (sum_BREAKF >ZERO)); 390 if (!breakFlag) { 391 breakFlag=FALSE; 392 /*** 393 ***** X_PRES and R_PRES are moved to memory: 394 ***/ 395 #pragma omp master 396 { 397 Factor=0.; 398 for (i=0;i0;--i) { 408 BREAKF[i]=BREAKF[i-1]; 409 R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1]; 410 R_PRES[i]=R_PRES[i-1]; 411 X_PRES[i]=X_PRES[i-1]; 412 P_PRES[i]=P_PRES[i-1]; 413 } 414 R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES; 415 R_PRES[0]=save_R_PRES; 416 X_PRES[0]=save_XPRES; 417 P_PRES[0]=save_P_PRES; 418 419 if (ABS(Factor)<=ZERO) { 420 Factor=1.; 421 BREAKF[0]=ZERO; 422 } else { 423 Factor=1./Factor; 424 BREAKF[0]=ONE; 425 } 426 for (i=0;i6) { 511 #pragma omp for private(z) schedule(static) 512 for (z=0;z0.) { 537 /*** 538 ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}: 539 ***/ 540 #pragma omp master 541 { 542 SC1=ZERO; 543 SC2=ZERO; 544 L2_R=ZERO; 545 } 546 #pragma omp barrier 547 #pragma omp for private(z,diff) reduction(+:SC1,SC2) schedule(static) 548 for (z=0;z0) restartFlag=(num_iter_restart >= restart); 563 } else { 564 convergeFlag=FALSE; 565 restartFlag=FALSE; 566 } 567 maxIterFlag = (num_iter >= maxit); 568 } 569 } 570 /* end of iteration */ 571 #pragma omp master 572 { 573 Norm_of_residual_global=norm_of_residual; 574 Num_iter_global=num_iter; 575 if (maxIterFlag) { 576 Status = SOLVER_MAXITER_REACHED; 577 } else if (breakFlag) { 578 Status = SOLVER_BREAKDOWN; 579 } 580 } 581 } /* end of parallel region */ 582 for (i=0;i

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision