# Diff of /branches/doubleplusgood/paso/src/TFQMR.cpp

revision 1711 by artak, Tue Aug 19 03:37:25 2008 UTC revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1
1
2  /*******************************************************  /*****************************************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2012 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * http://www.uq.edu.au
6   *  *
7   *                http://esscc.uq.edu.au  * Primary Business: Queensland, Australia
11   *  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12   *******************************************************/  * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15
16
17  /* TFQMR iterations */  /* TFQMR iterations */
18
# Line 22  Line 25
25  #include <omp.h>  #include <omp.h>
26  #endif  #endif
27
28  #ifdef PASO_MPI  #ifdef ESYS_MPI
29  #include <mpi.h>  #include <mpi.h>
30  #endif  #endif
31
# Line 40  Line 43
43  *  =========  *  =========
44  *  *
45  *  r       (input) DOUBLE PRECISION array, dimension N.  *  r       (input) DOUBLE PRECISION array, dimension N.
46  *          On entry, residual of inital guess x  *
47  *  *
48  *  x       (input/output) DOUBLE PRECISION array, dimension N.  *  x       (input/output) DOUBLE PRECISION array, dimension N.
49  *          On input, the initial guess.  *
50  *  *
51  *  ITER    (input/output) INT  *  ITER    (input/output) INT
52  *          On input, the maximum iterations to be performed.  *          On input, the maximum iterations to be performed.
# Line 52  Line 55
55  *  INFO    (output) INT  *  INFO    (output) INT
56  *  *
57  *          = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned.  *          = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned.
58  *          = SOLVEr_MAXITER_REACHED  *          = SOLVER_MAXITER_REACHED
59  *          = SOLVER_INPUT_ERROR Illegal parameter:  *          = SOLVER_INPUT_ERROR Illegal parameter:
60  *          = SOLVEr_BREAKDOWN: If parameters rHO or OMEGA become smaller  *          = SOLVER_BREAKDOWN: If parameters RHO or OMEGA become smaller
61  *          = SOLVER_MEMORY_ERROR : If parameters rHO or OMEGA become smaller  *          = SOLVER_MEMORY_ERROR : If parameters RHO or OMEGA become smaller
62  *  *
63  *  ==============================================================  *  ==============================================================
64  */  */
# Line 83  err_t Paso_Solver_TFQMR( Line 86  err_t Paso_Solver_TFQMR(
86    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
87    err_t status = SOLVER_NO_ERROR;    err_t status = SOLVER_NO_ERROR;
88    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
89    double  *u1=NULL, *u2=NULL, *y1=NULL, *y2=NULL, *d=NULL, *w=NULL, *v=NULL, *v_old=NULL, *tmp=NULL ;    double  *u1=NULL, *u2=NULL, *y1=NULL, *y2=NULL, *d=NULL, *w=NULL, *v=NULL, *v_old=NULL, *temp_vector=NULL,*res=NULL;
90
91    double eta,theta,tau,rho,beta,alpha,sigma,rhon,c;    double eta,theta,tau,rho,beta,alpha,sigma,rhon,c;
92
# Line 104  err_t Paso_Solver_TFQMR( Line 107  err_t Paso_Solver_TFQMR(
107    w=TMPMEMALLOC(n,double);    w=TMPMEMALLOC(n,double);
108    v=TMPMEMALLOC(n,double);    v=TMPMEMALLOC(n,double);
109    v_old=TMPMEMALLOC(n,double);    v_old=TMPMEMALLOC(n,double);
110        temp_vector=TMPMEMALLOC(n,double);
111        res=TMPMEMALLOC(n,double);
tmp=TMPMEMALLOC(n,double);

112
113   if (u1 ==NULL || u2== NULL || y1 == NULL || y2== NULL || d==NULL || w==NULL || v==NULL || v_old==NULL) {   if (u1 ==NULL || u2== NULL || y1 == NULL || y2== NULL || d==NULL || w==NULL || v==NULL || v_old==NULL) {
114       status=SOLVER_MEMORY_ERROR;       status=SOLVER_MEMORY_ERROR;
# Line 120  err_t Paso_Solver_TFQMR( Line 121  err_t Paso_Solver_TFQMR(
121      status=SOLVER_INPUT_ERROR;      status=SOLVER_INPUT_ERROR;
122    }    }
123
124      Paso_zeroes(n,x);
125
126    Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);    Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
127    Paso_Solver_solvePreconditioner(A,r,r);    Paso_SystemMatrix_solvePreconditioner(A,res,r);
128    Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);    Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
129
130    Performance_startMonitor(pp,PERFORMANCE_SOLVER);    Performance_startMonitor(pp,PERFORMANCE_SOLVER);
# Line 129  err_t Paso_Solver_TFQMR( Line 132  err_t Paso_Solver_TFQMR(
132    Paso_zeroes(n,u2);    Paso_zeroes(n,u2);
133    Paso_zeroes(n,y2);    Paso_zeroes(n,y2);
134
135    Paso_Copy(n,w,r);    Paso_Copy(n,w,res);
136    Paso_Copy(n,y1,r);    Paso_Copy(n,y1,res);
137
138    Paso_zeroes(n,d);    Paso_zeroes(n,d);
139
140    Performance_stopMonitor(pp,PERFORMANCE_SOLVER);    Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
141    Performance_startMonitor(pp,PERFORMANCE_MVM);    Performance_startMonitor(pp,PERFORMANCE_MVM);
142    Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, y1,ZERO,v);    Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, y1,PASO_ZERO,temp_vector);
143    Performance_stopMonitor(pp,PERFORMANCE_MVM);    Performance_stopMonitor(pp,PERFORMANCE_MVM);
144    Performance_startMonitor(pp,PERFORMANCE_SOLVER);    Performance_startMonitor(pp,PERFORMANCE_SOLVER);
145
146    Performance_stopMonitor(pp,PERFORMANCE_SOLVER);    Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
147    Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);    Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
148    Paso_Solver_solvePreconditioner(A,v,v);    Paso_SystemMatrix_solvePreconditioner(A,v,temp_vector);
149    Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);    Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
150    Performance_startMonitor(pp,PERFORMANCE_SOLVER);    Performance_startMonitor(pp,PERFORMANCE_SOLVER);
151
# Line 151  err_t Paso_Solver_TFQMR( Line 154  err_t Paso_Solver_TFQMR(
154    theta = 0.0;    theta = 0.0;
155    eta = 0.0;    eta = 0.0;
156
157    tau = Paso_l2(n,r,A->mpi_info);    tau = Paso_l2(n,res,A->mpi_info);
158
159    rho = tau * tau;    rho = tau * tau;
160
# Line 161  err_t Paso_Solver_TFQMR( Line 164  err_t Paso_Solver_TFQMR(
164    {    {
165
166
167       sigma=Paso_InnerProduct(n,r,v,A->mpi_info);       sigma=Paso_InnerProduct(n,res,v,A->mpi_info);
168
169       if (! (breakFlag = (ABS(sigma) == 0.))) {       if (! (breakFlag = (ABS(sigma) == 0.))) {
170
# Line 175  err_t Paso_Solver_TFQMR( Line 178  err_t Paso_Solver_TFQMR(
178
179            Performance_stopMonitor(pp,PERFORMANCE_SOLVER);            Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
180            Performance_startMonitor(pp,PERFORMANCE_MVM);            Performance_startMonitor(pp,PERFORMANCE_MVM);
181            Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, y2,ZERO,u2);            Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, y2,PASO_ZERO,temp_vector);
182            Performance_stopMonitor(pp,PERFORMANCE_MVM);            Performance_stopMonitor(pp,PERFORMANCE_MVM);
183            Performance_startMonitor(pp,PERFORMANCE_SOLVER);            Performance_startMonitor(pp,PERFORMANCE_SOLVER);
184
185            Performance_stopMonitor(pp,PERFORMANCE_SOLVER);            Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
186            Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);            Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
187            Paso_Solver_solvePreconditioner(A,u2,u2);        Paso_SystemMatrix_solvePreconditioner(A,u2,temp_vector);
188            Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);            Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
189            Performance_startMonitor(pp,PERFORMANCE_SOLVER);            Performance_startMonitor(pp,PERFORMANCE_SOLVER);
190           }           }
# Line 204  err_t Paso_Solver_TFQMR( Line 207  err_t Paso_Solver_TFQMR(
207
208       breakFlag = (ABS(rho) == 0);       breakFlag = (ABS(rho) == 0);
209
210       rhon = Paso_InnerProduct(n,r,w,A->mpi_info);       rhon = Paso_InnerProduct(n,res,w,A->mpi_info);
211       beta = rhon / rho;       beta = rhon / rho;
212       rho = rhon;       rho = rhon;
213
# Line 213  err_t Paso_Solver_TFQMR( Line 216  err_t Paso_Solver_TFQMR(
216       Performance_stopMonitor(pp,PERFORMANCE_SOLVER);       Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
217
218       Performance_startMonitor(pp,PERFORMANCE_MVM);       Performance_startMonitor(pp,PERFORMANCE_MVM);
219       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, y1,ZERO,u1);       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, y1,PASO_ZERO,temp_vector);
220       Performance_stopMonitor(pp,PERFORMANCE_MVM);       Performance_stopMonitor(pp,PERFORMANCE_MVM);
221
222       Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);       Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
223       Paso_Solver_solvePreconditioner(A,u1,u1);       Paso_SystemMatrix_solvePreconditioner(A,u1,temp_vector);
224       Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);       Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
225
226       Performance_startMonitor(pp,PERFORMANCE_SOLVER);       Performance_startMonitor(pp,PERFORMANCE_SOLVER);
# Line 250  err_t Paso_Solver_TFQMR( Line 253  err_t Paso_Solver_TFQMR(
253      TMPMEMFREE(w);      TMPMEMFREE(w);
254      TMPMEMFREE(v);      TMPMEMFREE(v);
255      TMPMEMFREE(v_old);      TMPMEMFREE(v_old);
256          TMPMEMFREE(temp_vector);
257        TMPMEMFREE(res);
258      *iter=num_iter;      *iter=num_iter;
259      *tolerance=norm_of_residual;      *tolerance=norm_of_residual;
260

Legend:
 Removed from v.1711 changed lines Added in v.3981