paso/src/Solver.c

/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* 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"
#include "esysUtils/blocktimer.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;
   double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b;
   double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local;
   double norm_max_of_residual_local,norm_max_of_residual;
   double last_norm_max_of_residual,*scaling;
#ifdef PASO_MPI
   double loc_norm;
#endif
   dim_t i,totIter,cntIter,method;
   bool_t finalizeIteration;
   err_t errorCode=NO_ERROR;
   dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
   dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
   double blocktimer_precond, blocktimer_start = blocktimer_time();
 
     Paso_resetError();
     tolerance=options->tolerance;
     if (tolerance < 100.* EPSILON) {
       Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small.");
     }
     if (tolerance >1.) {
       Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance mut be less than one.");
     }
     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric,A->mpi_info);
     /* check matrix type */
     if ((A->type & MATRIX_FORMAT_CSC) || (A->type & MATRIX_FORMAT_OFFSET1) || (A->type & MATRIX_FORMAT_SYM) ) {
       Paso_setError(TYPE_ERROR,"Paso_Solver: 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,"Paso_Solver: Iterative solver requires row and column block sizes to be equal.");
     }
     if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
        Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires a square matrix.");
        return;
     }
     /* this for testing only */
     if (method==PASO_NONLINEAR_GMRES) {
        Paso_Function* F=NULL;
        F=Paso_Function_LinearSystem_alloc(A,b,options);
        Paso_Solver_solvePreconditioner(A,x,b);
        errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
        if (errorCode!=NO_ERROR) {
           Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
        }
        Paso_Function_LinearSystem_free(F);
        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 < numEqua ; ++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);
                }
           }
           /* TODO: use one call */
           #ifdef PASO_MPI
           {
               loc_norm = norm2_of_b;
               MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
               loc_norm = norm_max_of_b;
               MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
           }
           #endif
           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 < numSol; i++) x[i]=0.;
            if (options->verbose) printf("right hand side is identical zero.\n");
         } else {
            if (options->verbose) {
               printf("Paso_Solver: l2/lmax-norm of right hand side is  %e/%e.\n",norm2_of_b,norm_max_of_b);
               printf("Paso_Solver: l2/lmax-stopping criterion is  %e/%e.\n",norm2_of_b*tolerance,norm_max_of_b*tolerance);
               switch (method) {
               case PASO_BICGSTAB:
                  printf("Paso_Solver: Iterative method is BiCGStab.\n");
                  break;
               case PASO_PCG:
                  printf("Paso_Solver: Iterative method is PCG.\n");
                  break;
               case PASO_TFQMR:
                  printf("Paso_Solver: Iterative method is TFQMR.\n");
                  break;
               case PASO_MINRES:
                  printf("Paso_Solver: Iterative method is MINRES.\n");
                  break;
               case PASO_PRES20:
                  printf("Paso_Solver: Iterative method is PRES20.\n");
                  break;
               case PASO_GMRES:
                  if (options->restart>0) {
                     printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
                  } else {
                     printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);
                  }
                  break;
               }
            }
            /* construct the preconditioner */
            blocktimer_precond = blocktimer_time();
            Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
            Paso_Solver_setPreconditioner(A,options);
            Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
            blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
            if (! Paso_noError()) return;
              time_iter=Paso_timer();
              /* get an initial guess by evaluating the preconditioner */
              Paso_Solver_solvePreconditioner(A,x,b);
              /* start the iteration process :*/
              r=TMPMEMALLOC(numEqua,double);
              Paso_checkPtr(r);
              if (Paso_noError()) {
                 totIter = 1;
                 finalizeIteration = FALSE;
                 last_norm2_of_residual=norm2_of_b;
                 last_norm_max_of_residual=norm_max_of_b;
                 /* Loop */
                 while (! finalizeIteration) {
                    cntIter = options->iter_max - totIter;
                    finalizeIteration = TRUE;
                    /*     Set initial residual. */
                    norm2_of_residual = 0;
                    norm_max_of_residual = 0;
                    #pragma omp parallel for private(i) schedule(static)
                    for (i = 0; i < numEqua; i++) r[i]=b[i];
                    Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
                    #pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
                    {
                       norm2_of_residual_local = 0;
                       norm_max_of_residual_local = 0;
                       #pragma omp for private(i) schedule(static)
                       for (i = 0; i < numEqua; 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);
                       }
                    }
                    /* TODO: use one call */
                    #ifdef PASO_MPI
                        loc_norm = norm2_of_residual;
                        MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
                        loc_norm = norm_max_of_residual;
                        MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
                    #endif
                    norm2_of_residual =sqrt(norm2_of_residual);
             
                  if (options->verbose) printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is  %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
                  if (totIter>1 && 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, "Paso_Solver: 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_TFQMR:
                           errorCode = Paso_Solver_TFQMR(A, r, x, &cntIter, &tol, pp);
                           break;
                        case PASO_MINRES:
                           errorCode = Paso_Solver_MINRES(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,"Paso_Solver: maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
                           if (options->verbose) printf("Paso_Solver: Maximum number of iterations reached.!\n");
                        } else if (errorCode == SOLVER_INPUT_ERROR ) {
                           Paso_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");
                           if (options->verbose) printf("Paso_Solver: Internal error!\n");
                        } else if ( errorCode == SOLVER_BREAKDOWN ) {
                           if (cntIter <= 1) {
                              Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");
                              if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");
                           } else {
                              if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
                              finalizeIteration = FALSE;
                           }
                      } else {
                         if (options->verbose) printf(". convergence! \n");
                      }
                   }
                 } /* while */
              }
              MEMFREE(r);
              time_iter=Paso_timer()-time_iter;
              if (options->verbose)  {
                 printf("\ntiming: Paso_Solver:  %.4e sec\n",time_iter);
                 if (totIter>1) {
                    if(totIter==options->iter_max) {
                        printf("timing: Total MAX steps, time per iteration step: %.4e sec\n",time_iter/totIter);
                    } else {
                        printf("timing: Total %d steps, time per iteration step: %.4e sec\n",totIter,time_iter/totIter);    
                    }
                 }
                 else {
                    printf("timing: Total 1 step, time per iteration step: %.4e sec\n",time_iter);
                    }   
              }
           }
        }
      }
   }
   Performance_stopMonitor(pp,PERFORMANCE_ALL);
   blocktimer_increment("Paso_Solver()", blocktimer_start);
}
1	/*******************************************************
2	*
3	* Copyright (c) 2003-2008 by University of Queensland
4	* Earth Systems Science Computational Center (ESSCC)
5	* http://www.uq.edu.au/esscc
6	*
7	* Primary Business: Queensland, Australia
8	* Licensed under the Open Software License version 3.0
9	* http://www.opensource.org/licenses/osl-3.0.php
10	*
11	*******************************************************/
12
13
14	/**************************************************************/
15
16	/* Paso: SystemMatrix: controls iterative linear system solvers */
17
18	/**************************************************************/
19
20	/* Copyrights by ACcESS Australia 2003/04 */
21	/* Author: gross@access.edu.au */
22
23	/**************************************************************/
24
25	#include "Paso.h"
26	#include "SystemMatrix.h"
27	#include "Solver.h"
28	#include "esysUtils/blocktimer.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,
41	Paso_Options* options,Paso_Performance* pp) {
42
43	double norm2_of_b,tol,tolerance,time_iter;
44	double *r=NULL,norm2_of_residual,last_norm2_of_residual,norm_max_of_b;
45	double norm2_of_b_local,norm_max_of_b_local,norm2_of_residual_local;
46	double norm_max_of_residual_local,norm_max_of_residual;
47	double last_norm_max_of_residual,*scaling;
48	#ifdef PASO_MPI
49	double loc_norm;
50	#endif
51	dim_t i,totIter,cntIter,method;
52	bool_t finalizeIteration;
53	err_t errorCode=NO_ERROR;
54	dim_t numSol = Paso_SystemMatrix_getTotalNumCols(A);
55	dim_t numEqua = Paso_SystemMatrix_getTotalNumRows(A);
56	double blocktimer_precond, blocktimer_start = blocktimer_time();
57
58	Paso_resetError();
59	tolerance=options->tolerance;
60	if (tolerance < 100.* EPSILON) {
61	Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance is too small.");
62	}
63	if (tolerance >1.) {
64	Paso_setError(VALUE_ERROR,"Paso_Solver: Tolerance mut be less than one.");
65	}
66	method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric,A->mpi_info);
67	/* check matrix type */
68	if ((A->type & MATRIX_FORMAT_CSC) \|\| (A->type & MATRIX_FORMAT_OFFSET1) \|\| (A->type & MATRIX_FORMAT_SYM) ) {
69	Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires CSR format with unsymmetric storage scheme and index offset 0.");
70	}
71	if (A->col_block_size != A->row_block_size) {
72	Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires row and column block sizes to be equal.");
73	}
74	if (Paso_SystemMatrix_getGlobalNumCols(A) != Paso_SystemMatrix_getGlobalNumRows(A)) {
75	Paso_setError(TYPE_ERROR,"Paso_Solver: Iterative solver requires a square matrix.");
76	return;
77	}
78	/* this for testing only */
79	if (method==PASO_NONLINEAR_GMRES) {
80	Paso_Function* F=NULL;
81	F=Paso_Function_LinearSystem_alloc(A,b,options);
82	Paso_Solver_solvePreconditioner(A,x,b);
83	errorCode=Paso_Solver_NewtonGMRES(F,x,options,pp);
84	if (errorCode!=NO_ERROR) {
85	Paso_setError(SYSTEM_ERROR,"Paso_Solver_NewtonGMRES: an error has occured.");
86	}
87	Paso_Function_LinearSystem_free(F);
88	return;
89	}
90	/* ========================= */
91	Performance_startMonitor(pp,PERFORMANCE_ALL);
92	if (Paso_noError()) {
93	/* get normalization */
94	scaling=Paso_SystemMatrix_borrowNormalization(A);
95	if (Paso_noError()) {
96	/* get the norm of the right hand side */
97	norm2_of_b=0.;
98	norm_max_of_b=0.;
99	#pragma omp parallel private(norm2_of_b_local,norm_max_of_b_local)
100	{
101	norm2_of_b_local=0.;
102	norm_max_of_b_local=0.;
103	#pragma omp for private(i) schedule(static)
104	for (i = 0; i < numEqua ; ++i) {
105	norm2_of_b_local += b[i] * b[i];
106	norm_max_of_b_local = MAX(ABS(scaling[i]*b[i]),norm_max_of_b_local);
107	}
108	#pragma omp critical
109	{
110	norm2_of_b += norm2_of_b_local;
111	norm_max_of_b = MAX(norm_max_of_b_local,norm_max_of_b);
112	}
113	}
114	/* TODO: use one call */
115	#ifdef PASO_MPI
116	{
117	loc_norm = norm2_of_b;
118	MPI_Allreduce(&loc_norm,&norm2_of_b, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
119	loc_norm = norm_max_of_b;
120	MPI_Allreduce(&loc_norm,&norm_max_of_b, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
121	}
122	#endif
123	norm2_of_b=sqrt(norm2_of_b);
124
125	/* if norm2_of_b==0 we are ready: x=0 */
126	if (norm2_of_b <=0.) {
127	#pragma omp parallel for private(i) schedule(static)
128	for (i = 0; i < numSol; i++) x[i]=0.;
129	if (options->verbose) printf("right hand side is identical zero.\n");
130	} else {
131	if (options->verbose) {
132	printf("Paso_Solver: l2/lmax-norm of right hand side is %e/%e.\n",norm2_of_b,norm_max_of_b);
133	printf("Paso_Solver: l2/lmax-stopping criterion is %e/%e.\n",norm2_of_btolerance,norm_max_of_btolerance);
134	switch (method) {
135	case PASO_BICGSTAB:
136	printf("Paso_Solver: Iterative method is BiCGStab.\n");
137	break;
138	case PASO_PCG:
139	printf("Paso_Solver: Iterative method is PCG.\n");
140	break;
141	case PASO_TFQMR:
142	printf("Paso_Solver: Iterative method is TFQMR.\n");
143	break;
144	case PASO_MINRES:
145	printf("Paso_Solver: Iterative method is MINRES.\n");
146	break;
147	case PASO_PRES20:
148	printf("Paso_Solver: Iterative method is PRES20.\n");
149	break;
150	case PASO_GMRES:
151	if (options->restart>0) {
152	printf("Paso_Solver: Iterative method is GMRES(%d,%d)\n",options->truncation,options->restart);
153	} else {
154	printf("Paso_Solver: Iterative method is GMRES(%d)\n",options->truncation);
155	}
156	break;
157	}
158	}
159	/* construct the preconditioner */
160	blocktimer_precond = blocktimer_time();
161	Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
162	Paso_Solver_setPreconditioner(A,options);
163	Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER_INIT);
164	blocktimer_increment("Paso_Solver_setPreconditioner()", blocktimer_precond);
165	if (! Paso_noError()) return;
166	time_iter=Paso_timer();
167	/* get an initial guess by evaluating the preconditioner */
168	Paso_Solver_solvePreconditioner(A,x,b);
169	/* start the iteration process :*/
170	r=TMPMEMALLOC(numEqua,double);
171	Paso_checkPtr(r);
172	if (Paso_noError()) {
173	totIter = 1;
174	finalizeIteration = FALSE;
175	last_norm2_of_residual=norm2_of_b;
176	last_norm_max_of_residual=norm_max_of_b;
177	/* Loop */
178	while (! finalizeIteration) {
179	cntIter = options->iter_max - totIter;
180	finalizeIteration = TRUE;
181	/* Set initial residual. */
182	norm2_of_residual = 0;
183	norm_max_of_residual = 0;
184	#pragma omp parallel for private(i) schedule(static)
185	for (i = 0; i < numEqua; i++) r[i]=b[i];
186	Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), r);
187	#pragma omp parallel private(norm2_of_residual_local,norm_max_of_residual_local)
188	{
189	norm2_of_residual_local = 0;
190	norm_max_of_residual_local = 0;
191	#pragma omp for private(i) schedule(static)
192	for (i = 0; i < numEqua; i++) {
193	norm2_of_residual_local+= r[i] * r[i];
194	norm_max_of_residual_local=MAX(ABS(scaling[i]*r[i]),norm_max_of_residual_local);
195	}
196	#pragma omp critical
197	{
198	norm2_of_residual += norm2_of_residual_local;
199	norm_max_of_residual = MAX(norm_max_of_residual_local,norm_max_of_residual);
200	}
201	}
202	/* TODO: use one call */
203	#ifdef PASO_MPI
204	loc_norm = norm2_of_residual;
205	MPI_Allreduce(&loc_norm,&norm2_of_residual, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
206	loc_norm = norm_max_of_residual;
207	MPI_Allreduce(&loc_norm,&norm_max_of_residual, 1, MPI_DOUBLE, MPI_MAX, A->mpi_info->comm);
208	#endif
209	norm2_of_residual =sqrt(norm2_of_residual);
210
211	if (options->verbose) printf("Paso_Solver: Step %5d: l2/lmax-norm of residual is %e/%e",totIter,norm2_of_residual,norm_max_of_residual);
212	if (totIter>1 && norm2_of_residual>=last_norm2_of_residual && norm_max_of_residual>=last_norm_max_of_residual) {
213	if (options->verbose) printf(" divergence!\n");
214	Paso_setError(WARNING, "Paso_Solver: No improvement during iteration. Iterative solver gives up.");
215	} else {
216	/* */
217	if (norm2_of_residual>tolerancenorm2_of_b \|\| norm_max_of_residual>tolerancenorm_max_of_b ) {
218	tol=toleranceMIN(norm2_of_b,0.1norm2_of_residual/norm_max_of_residual*norm_max_of_b);
219	if (options->verbose) printf(" (new tolerance = %e).\n",tol);
220	last_norm2_of_residual=norm2_of_residual;
221	last_norm_max_of_residual=norm_max_of_residual;
222	/* call the solver */
223	switch (method) {
224	case PASO_BICGSTAB:
225	errorCode = Paso_Solver_BiCGStab(A, r, x, &cntIter, &tol, pp);
226	break;
227	case PASO_PCG:
228	errorCode = Paso_Solver_PCG(A, r, x, &cntIter, &tol, pp);
229	break;
230	case PASO_TFQMR:
231	errorCode = Paso_Solver_TFQMR(A, r, x, &cntIter, &tol, pp);
232	break;
233	case PASO_MINRES:
234	errorCode = Paso_Solver_MINRES(A, r, x, &cntIter, &tol, pp);
235	break;
236	case PASO_PRES20:
237	errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,5,20, pp);
238	break;
239	case PASO_GMRES:
240	errorCode = Paso_Solver_GMRES(A, r, x, &cntIter, &tol,options->truncation,options->restart, pp);
241	break;
242	}
243	totIter += cntIter;
244	/* error handling */
245	if (errorCode==NO_ERROR) {
246	finalizeIteration = FALSE;
247	} else if (errorCode==SOLVER_MAXITER_REACHED) {
248	Paso_setError(WARNING,"Paso_Solver: maximum number of iteration step reached.\nReturned solution does not fulfil stopping criterion.");
249	if (options->verbose) printf("Paso_Solver: Maximum number of iterations reached.!\n");
250	} else if (errorCode == SOLVER_INPUT_ERROR ) {
251	Paso_setError(SYSTEM_ERROR,"Paso_Solver: illegal dimension in iterative solver.");
252	if (options->verbose) printf("Paso_Solver: Internal error!\n");
253	} else if ( errorCode == SOLVER_BREAKDOWN ) {
254	if (cntIter <= 1) {
255	Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver: fatal break down in iterative solver.");
256	if (options->verbose) printf("Paso_Solver: Uncurable break down!\n");
257	} else {
258	if (options->verbose) printf("Paso_Solver: Breakdown at iter %d (residual = %e). Restarting ...\n", cntIter+totIter, tol);
259	finalizeIteration = FALSE;
260	}
261	} else {
262	if (options->verbose) printf(". convergence! \n");
263	}
264	}
265	} /* while */
266	}
267	MEMFREE(r);
268	time_iter=Paso_timer()-time_iter;
269	if (options->verbose) {
270	printf("\ntiming: Paso_Solver: %.4e sec\n",time_iter);
271	if (totIter>1) {
272	if(totIter==options->iter_max) {
273	printf("timing: Total MAX steps, time per iteration step: %.4e sec\n",time_iter/totIter);
274	} else {
275	printf("timing: Total %d steps, time per iteration step: %.4e sec\n",totIter,time_iter/totIter);
276	}
277	}
278	else {
279	printf("timing: Total 1 step, time per iteration step: %.4e sec\n",time_iter);
280	}
281	}
282	}
283	}
284	}
285	}
286	Performance_stopMonitor(pp,PERFORMANCE_ALL);
287	blocktimer_increment("Paso_Solver()", blocktimer_start);
288	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision