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

Revision 633 - (show annotations)
Thu Mar 23 05:37:00 2006 UTC (13 years, 11 months ago) by dhawcroft
Original Path: trunk/paso/src/Solvers/GMRES.c
File MIME type: text/plain
File size: 23862 byte(s)

 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[MAX(Length_of_recursion,0)+1],*R_PRES[MAX(Length_of_recursion,0)+1],*P_PRES[MAX(Length_of_recursion,0)+1]; 77 double P_PRES_dot_AP[MAX(Length_of_recursion,0)],R_PRES_dot_P_PRES[MAX(Length_of_recursion,0)+1],BREAKF[MAX(Length_of_recursion,0)+1],ALPHA[MAX(Length_of_recursion,0)]; 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; 83 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE; 84 err_t Status=SOLVER_NO_ERROR; 85 86 /* adapt original routine parameters */ 87 88 dim_t n=A->num_cols * A-> col_block_size; 89 dim_t Length_of_mem=MAX(Length_of_recursion,0)+1; 90 91 /* Test the input parameters. */ 92 if (restart>0) restart=MAX(Length_of_recursion,restart); 93 if (n < 0 || Length_of_recursion<=0) { 94 return SOLVER_INPUT_ERROR; 95 } 96 97 /* allocate memory: */ 98 99 AP=TMPMEMALLOC(n,double); 100 if (AP==NULL) Status=SOLVER_MEMORY_ERROR; 101 for (i=0;i6) { 320 #pragma omp master 321 { 322 R_PRES_dot_P=ZERO; 323 P_PRES_dot_AP0=ZERO; 324 P_PRES_dot_AP1=ZERO; 325 P_PRES_dot_AP2=ZERO; 326 P_PRES_dot_AP3=ZERO; 327 P_PRES_dot_AP4=ZERO; 328 P_PRES_dot_AP5=ZERO; 329 P_PRES_dot_AP6=ZERO; 330 } 331 #pragma omp barrier 332 #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) 333 for (z=0;z ZERO) && (sum_BREAKF >ZERO)); 378 if (!breakFlag) { 379 breakFlag=FALSE; 380 /*** 381 ***** X_PRES and R_PRES are moved to memory: 382 ***/ 383 #pragma omp master 384 { 385 Factor=0.; 386 for (i=0;i0;--i) { 396 BREAKF[i]=BREAKF[i-1]; 397 R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1]; 398 R_PRES[i]=R_PRES[i-1]; 399 X_PRES[i]=X_PRES[i-1]; 400 P_PRES[i]=P_PRES[i-1]; 401 } 402 R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES; 403 R_PRES[0]=save_R_PRES; 404 X_PRES[0]=save_XPRES; 405 P_PRES[0]=save_P_PRES; 406 407 if (ABS(Factor)<=ZERO) { 408 Factor=1.; 409 BREAKF[0]=ZERO; 410 } else { 411 Factor=1./Factor; 412 BREAKF[0]=ONE; 413 } 414 for (i=0;i6) { 499 #pragma omp for private(z) schedule(static) 500 for (z=0;z0.) { 525 /*** 526 ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}: 527 ***/ 528 #pragma omp master 529 { 530 SC1=ZERO; 531 SC2=ZERO; 532 L2_R=ZERO; 533 } 534 #pragma omp barrier 535 #pragma omp for private(z,diff) reduction(+:SC1,SC2) schedule(static) 536 for (z=0;z0) restartFlag=(num_iter_restart >= restart); 551 } else { 552 convergeFlag=FALSE; 553 restartFlag=FALSE; 554 } 555 maxIterFlag = (num_iter >= maxit); 556 } 557 } 558 /* end of iteration */ 559 #pragma omp master 560 { 561 Norm_of_residual_global=norm_of_residual; 562 Num_iter_global=num_iter; 563 if (maxIterFlag) { 564 Status = SOLVER_MAXITER_REACHED; 565 } else if (breakFlag) { 566 Status = SOLVER_BREAKDOWN; 567 } 568 } 569 } /* end of parallel region */ 570 TMPMEMFREE(AP); 571 for (i=0;i

## Properties

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