src/Solvers/Solver.c

/* $Id$ */


/*
********************************************************************************
*               Copyright   2006 by ACcESS MNRF                                *
*                                                                              * 
*                 http://www.access.edu.au                                     *
*           Primary Business: Queensland, Australia                            *
*     Licensed under the Open Software License version 3.0             *
*        http://www.opensource.org/licenses/osl-3.0.php                        *
********************************************************************************
*/

/**************************************************************/

/* Paso: SystemMatrix: controls iterative linear system solvers  */ 

/**************************************************************/

/* Copyrights by ACcESS Australia 2003/04 */
/* Author: gross@access.edu.au */

/**************************************************************/

#include "../Paso.h"
#include "../SystemMatrix.h"
#include "Solver.h"

/***********************************************************************************/

/*  free space */

void Paso_Solver_free(Paso_SystemMatrix* A) {
    Paso_Preconditioner_free(A->solver);
    A->solver=NULL;
}
/*  call the iterative solver: */

void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) {
    double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,
           norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual,
           last_norm_max_of_residual,*scaling;
     dim_t i,totIter,cntIter,method;
     bool_t finalizeIteration;
     err_t errorCode;
     dim_t n_col = A->num_cols * A-> col_block_size;
     dim_t n_row = A->num_rows * A-> row_block_size;
 
     tolerance=MAX(options->tolerance,EPSILON);
     Paso_resetError();
     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
     /* check matrix type */
     if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
       Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
     }
     if (A->col_block_size != A->row_block_size) {
        Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");
     }
     if (A->num_cols != A->num_rows) {
        Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
        return;
     }
     Performance_startMonitor(pp,PERFORMANCE_ALL);
     if (Paso_noError()) {
        /* get normalization */
        scaling=Paso_SystemMatrix_borrowNormalization(A);
        if (Paso_noError()) {
           /* get the norm of the right hand side */
           norm2_of_b=0.;
           norm_max_of_b=0.;
           #pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local) 
           { 
                norm2_of_b_local=0.;
                norm_max_of_b_local=0.;
                #pragma omp for private(i) schedule(static)
                for (i = 0; i < n_row ; ++i) {
                      norm2_of_b_local += b[i] * b[i];
                      norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
                }
                #pragma omp critical
                {
                   norm2_of_b += norm2_of_b_local;
                   norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
                }
           }
           norm2_of_b=sqrt(norm2_of_b);
    
           /* if norm2_of_b==0 we are ready: x=0 */
           if (norm2_of_b <=0.) {
               #pragma omp parallel for private(i) schedule(static)
               for (i = 0; i < n_col; i++) x[i]=0.;
               if (options->verbose) printf("right hand side is identical zero.\n");
           } else {
              if (options->verbose) {
                  printf("l2/lmax-norm of right hand side is  %e/%e.\n",norm2_of_b,norm_max_of_b);
                  printf("l2/lmax-stopping criterion is  %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
                  switch (method) {
                     case PASO_BICGSTAB:
                         printf("Iterative method is BiCGStab.\n");
                         break;
                     case PASO_PCG:
                         printf("Iterative method is PCG.\n");
                         break;
                     case PASO_PRES20:
                         printf("Iterative method is PRES20.\n");
                         break;
                     case PASO_GMRES:
                         if (options->restart>0) {
                            printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
                         } else {
                            printf("Iterative method is GMRES(%d)\n",options->truncation);
                         }
                         break;
                  }
              }
              /* construct the preconditioner */
              
              Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
              Paso_Solver_setPreconditioner(A,options);
              Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
              if (! Paso_noError()) return;
          
              /* get an initial guess by evaluating the preconditioner */
              time_iter=Paso_timer();
              #pragma omp parallel
              Paso_Solver_solvePreconditioner(A,x,b);
              if (! Paso_noError()) return;
              /* start the iteration process :*/
              r=TMPMEMALLOC(n_row,double);
              if (Paso_checkPtr(r)) return;
              totIter = 0;
              finalizeIteration = FALSE;
              last_norm2_of_residual=norm2_of_b;
              last_norm_max_of_residual=norm_max_of_b;
              while (! finalizeIteration) {
                 cntIter = options->iter_max - totIter;
                 finalizeIteration = TRUE;
                 /*     Set initial residual. */
                 norm2_of_residual = 0;
                 norm_max_of_residual = 0;
                 #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
                 {
                    #pragma omp for private(i) schedule(static)
                    for (i = 0; i < n_row; i++) r[i]=b[i];
          
                    Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
          
                    norm2_of_residual_local = 0;
                    norm_max_of_residual_local = 0;
                    #pragma omp for private(i) schedule(static)
                    for (i = 0; i < n_row; i++) {
                           norm2_of_residual_local+= r[i] * r[i];
                           norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
                    }
                    #pragma omp critical
                    {
                       norm2_of_residual += norm2_of_residual_local;
                       norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
                    }
                 }
                 norm2_of_residual =sqrt(norm2_of_residual);
          
                 if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
                 if (totIter>0 && norm2_of_residual>=last_norm2_of_residual &&  norm_max_of_residual>=last_norm_max_of_residual) {
                      if (options->verbose) printf(" divergence!\n");
                      Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
                 } else {
                    /* */
                    if (norm2_of_residual>tolerance*norm2_of_b || norm_max_of_residual>tolerance*norm_max_of_b ) {
                        tol=tolerance*MIN(norm2_of_b,0.1*norm2_of_residual/norm_max_of_residual*norm_max_of_b);
                       if (options->verbose) printf(" (new tolerance = %e).\n",tol);
                       last_norm2_of_residual=norm2_of_residual;
                       last_norm_max_of_residual=norm_max_of_residual;
                       /* call the solver */
                       switch (method) {
                              case PASO_BICGSTAB:
                                  errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol,pp); 
                                  break;
                              case PASO_PCG:
                                  errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol,pp); 
                                  break;
                              case PASO_PRES20:
                                  errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20,pp); 
                                  break;
                              case PASO_GMRES:
                                  errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart,pp); 
                                  break;
                        }
                        totIter += cntIter;
                   
                        /* error handling  */
                        if (errorCode==NO_ERROR) {
                              finalizeIteration = FALSE;
                         } else if (errorCode==SOLVER_MAXITER_REACHED) {
                               Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
                               if (options->verbose) printf("Maximum number of iterations reached.!\n");
                        } else if (errorCode == SOLVER_INPUT_ERROR ) {
                               Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
                               if (options->verbose) printf("Internal error!\n");
                        } else if ( errorCode == SOLVER_BREAKDOWN ) {
                         if (cntIter <= 1) {
                                 Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
                                 if (options->verbose) printf("Uncurable break down!\n");
                            } else {
                             if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
                                finalizeIteration = FALSE;
                            }
                        } else if (errorCode == SOLVER_MEMORY_ERROR) {
                          Paso_setError(MEMORY_ERROR,"memory allocation failed.");
                             if (options->verbose) printf("Memory allocation failed!\n");
                        } else if (errorCode !=SOLVER_NO_ERROR ) {
                          Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
                             if (options->verbose) printf("Unidentified error!\n");
                        }
                   } else {
                      if (options->verbose) printf(". convergence! \n");
                   }
                }
              } /* while */
              MEMFREE(r);
              time_iter=Paso_timer()-time_iter;
              if (options->verbose)  {
                 printf("timing: solver: %.4e sec\n",time_iter);
                 if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
              }
           }
        }
      }
      Performance_stopMonitor(pp,PERFORMANCE_ALL);
}
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
17	/* Paso: SystemMatrix: controls iterative linear system solvers */
18
19	/**************************************************************/
20
21	/* Copyrights by ACcESS Australia 2003/04 */
22	/* Author: gross@access.edu.au */
23
24	/**************************************************************/
25
26	#include "../Paso.h"
27	#include "../SystemMatrix.h"
28	#include "Solver.h"
29
30	/***********************************************************************************/
31
32	/* free space */
33
34	void Paso_Solver_free(Paso_SystemMatrix* A) {
35	Paso_Preconditioner_free(A->solver);
36	A->solver=NULL;
37	}
38	/* call the iterative solver: */
39
40	void Paso_Solver(Paso_SystemMatrix* A,double* x,double* b,Paso_Options* options,Paso_Performance* pp) {
41	double norm2_of_b,tol,tolerance,time_iter,*r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b,
42	norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local,norm_max_of_residual_local,norm_max_of_residual,
43	last_norm_max_of_residual,*scaling;
44	dim_t i,totIter,cntIter,method;
45	bool_t finalizeIteration;
46	err_t errorCode;
47	dim_t n_col = A->num_cols * A-> col_block_size;
48	dim_t n_row = A->num_rows * A-> row_block_size;
49
50	tolerance=MAX(options->tolerance,EPSILON);
51	Paso_resetError();
52	method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
53	/* check matrix type */
54	if ((A->type & MATRIX_FORMAT_CSC) \|\| (A->type & MATRIX_FORMAT_OFFSET1) \|\| (A->type & MATRIX_FORMAT_SYM) ) {
55	Paso_setError(TYPE_ERROR,"Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
56	}
57	if (A->col_block_size != A->row_block_size) {
58	Paso_setError(TYPE_ERROR,"Iterative solver requires row and column block sizes to be equal.");
59	}
60	if (A->num_cols != A->num_rows) {
61	Paso_setError(TYPE_ERROR,"Iterative solver requires requires a square matrix.");
62	return;
63	}
64	Performance_startMonitor(pp,PERFORMANCE_ALL);
65	if (Paso_noError()) {
66	/* get normalization */
67	scaling=Paso_SystemMatrix_borrowNormalization(A);
68	if (Paso_noError()) {
69	/* get the norm of the right hand side */
70	norm2_of_b=0.;
71	norm_max_of_b=0.;
72	#pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
73	{
74	norm2_of_b_local=0.;
75	norm_max_of_b_local=0.;
76	#pragma omp for private(i) schedule(static)
77	for (i = 0; i < n_row ; ++i) {
78	norm2_of_b_local += b[i] * b[i];
79	norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
80	}
81	#pragma omp critical
82	{
83	norm2_of_b += norm2_of_b_local;
84	norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
85	}
86	}
87	norm2_of_b=sqrt(norm2_of_b);
88
89	/* if norm2_of_b==0 we are ready: x=0 */
90	if (norm2_of_b <=0.) {
91	#pragma omp parallel for private(i) schedule(static)
92	for (i = 0; i < n_col; i++) x[i]=0.;
93	if (options->verbose) printf("right hand side is identical zero.\n");
94	} else {
95	if (options->verbose) {
96	printf("l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
97	printf("l2/lmax-stopping criterion is %e/%e.\n",norm2_of_btolerance,norm_max_of_btolerance);
98	switch (method) {
99	case PASO_BICGSTAB:
100	printf("Iterative method is BiCGStab.\n");
101	break;
102	case PASO_PCG:
103	printf("Iterative method is PCG.\n");
104	break;
105	case PASO_PRES20:
106	printf("Iterative method is PRES20.\n");
107	break;
108	case PASO_GMRES:
109	if (options->restart>0) {
110	printf("Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
111	} else {
112	printf("Iterative method is GMRES(%d)\n",options->truncation);
113	}
114	break;
115	}
116	}
117	/* construct the preconditioner */
118
119	Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
120	Paso_Solver_setPreconditioner(A,options);
121	Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
122	if (! Paso_noError()) return;
123
124	/* get an initial guess by evaluating the preconditioner */
125	time_iter=Paso_timer();
126	#pragma omp parallel
127	Paso_Solver_solvePreconditioner(A,x,b);
128	if (! Paso_noError()) return;
129	/* start the iteration process :*/
130	r=TMPMEMALLOC(n_row,double);
131	if (Paso_checkPtr(r)) return;
132	totIter = 0;
133	finalizeIteration = FALSE;
134	last_norm2_of_residual=norm2_of_b;
135	last_norm_max_of_residual=norm_max_of_b;
136	while (! finalizeIteration) {
137	cntIter = options->iter_max - totIter;
138	finalizeIteration = TRUE;
139	/* Set initial residual. */
140	norm2_of_residual = 0;
141	norm_max_of_residual = 0;
142	#pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
143	{
144	#pragma omp for private(i) schedule(static)
145	for (i = 0; i < n_row; i++) r[i]=b[i];
146
147	Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
148
149	norm2_of_residual_local = 0;
150	norm_max_of_residual_local = 0;
151	#pragma omp for private(i) schedule(static)
152	for (i = 0; i < n_row; i++) {
153	norm2_of_residual_local+= r[i] * r[i];
154	norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
155	}
156	#pragma omp critical
157	{
158	norm2_of_residual += norm2_of_residual_local;
159	norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
160	}
161	}
162	norm2_of_residual =sqrt(norm2_of_residual);
163
164	if (options->verbose) printf("Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
165	if (totIter>0 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
166	if (options->verbose) printf(" divergence!\n");
167	Paso_setError(WARNING, "No improvement during iteration. Iterative solver gives up.");
168	} else {
169	/* */
170	if (norm2_of_residual>tolerancenorm2_of_b \|\| norm_max_of_residual>tolerancenorm_max_of_b ) {
171	tol=toleranceMIN(norm2_of_b,0.1norm2_of_residual/norm_max_of_residual*norm_max_of_b);
172	if (options->verbose) printf(" (new tolerance = %e).\n",tol);
173	last_norm2_of_residual=norm2_of_residual;
174	last_norm_max_of_residual=norm_max_of_residual;
175	/* call the solver */
176	switch (method) {
177	case PASO_BICGSTAB:
178	errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol,pp);
179	break;
180	case PASO_PCG:
181	errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol,pp);
182	break;
183	case PASO_PRES20:
184	errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20,pp);
185	break;
186	case PASO_GMRES:
187	errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart,pp);
188	break;
189	}
190	totIter += cntIter;
191
192	/* error handling */
193	if (errorCode==NO_ERROR) {
194	finalizeIteration = FALSE;
195	} else if (errorCode==SOLVER_MAXITER_REACHED) {
196	Paso_setError(WARNING,"maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
197	if (options->verbose) printf("Maximum number of iterations reached.!\n");
198	} else if (errorCode == SOLVER_INPUT_ERROR ) {
199	Paso_setError(SYSTEM_ERROR,"illegal dimension in iterative solver.");
200	if (options->verbose) printf("Internal error!\n");
201	} else if ( errorCode == SOLVER_BREAKDOWN ) {
202	if (cntIter <= 1) {
203	Paso_setError(ZERO_DIVISION_ERROR, "fatal break down in iterative solver.");
204	if (options->verbose) printf("Uncurable break down!\n");
205	} else {
206	if (options->verbose) printf("Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
207	finalizeIteration = FALSE;
208	}
209	} else if (errorCode == SOLVER_MEMORY_ERROR) {
210	Paso_setError(MEMORY_ERROR,"memory allocation failed.");
211	if (options->verbose) printf("Memory allocation failed!\n");
212	} else if (errorCode !=SOLVER_NO_ERROR ) {
213	Paso_setError(SYSTEM_ERROR,"unidentified error in iterative solver.");
214	if (options->verbose) printf("Unidentified error!\n");
215	}
216	} else {
217	if (options->verbose) printf(". convergence! \n");
218	}
219	}
220	} /* while */
221	MEMFREE(r);
222	time_iter=Paso_timer()-time_iter;
223	if (options->verbose) {
224	printf("timing: solver: %.4e sec\n",time_iter);
225	if (totIter>0) printf("timing: per iteration step: %.4e sec\n",time_iter/totIter);
226	}
227	}
228	}
229	}
230	Performance_stopMonitor(pp,PERFORMANCE_ALL);
231	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision