PaStiX Handbook  6.3.2
c_refine_grad.c
Go to the documentation of this file.
1 /**
2  *
3  * @file c_refine_grad.c
4  *
5  * PaStiX refinement functions implementations.
6  *
7  * @copyright 2015-2023 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8  * Univ. Bordeaux. All rights reserved.
9  *
10  * @version 6.3.2
11  * @author Mathieu Faverge
12  * @author Pierre Ramet
13  * @author Xavier Lacoste
14  * @author Theophile Terraz
15  * @author Gregoire Pichon
16  * @author Vincent Bridonneau
17  * @date 2023-07-21
18  * @generated from /builds/solverstack/pastix/refinement/z_refine_grad.c, normal z -> c, Wed Dec 13 12:09:46 2023
19  *
20  **/
21 #include "common.h"
22 #include "bcsc/bcsc.h"
23 #include "c_refine_functions.h"
24 
25 /**
26  *******************************************************************************
27  *
28  * @ingroup pastix_refine
29  *
30  * c_grad_smp - Refine the solution using conjugate gradian method.
31  *
32  *******************************************************************************
33  *
34  * @param[in] pastix_data
35  * The PaStiX data structure that describes the solver instance.
36  *
37  * @param[out] xp
38  * The solution vector.
39  *
40  * @param[in] bp
41  * The right hand side member (only one).
42  *
43  *******************************************************************************
44  *
45  * @return Number of iterations
46  *
47  *******************************************************************************/
49 c_grad_smp( pastix_data_t *pastix_data,
50  pastix_rhs_t xp,
51  pastix_rhs_t bp )
52 {
53  struct c_solver solver;
54  pastix_int_t n;
55  Clock refine_clk;
56  pastix_fixdbl_t t0 = 0;
57  pastix_fixdbl_t t3 = 0;
58  int itermax;
59  int nb_iter = 0;
60  int precond = 1;
63  pastix_complex32_t *gradr;
64  pastix_complex32_t *gradp;
65  pastix_complex32_t *gradz;
66  pastix_complex32_t *grad2;
67  pastix_complex32_t *sgrad = NULL;
68  float normb, normx, normr, alpha, beta;
69  float resid_b, eps;
70 
71  memset( &solver, 0, sizeof(struct c_solver) );
72  c_refine_init( &solver, pastix_data );
73 
74  if ( !(pastix_data->steps & STEP_NUMFACT) ) {
75  precond = 0;
76  }
77 
78  n = pastix_data->bcsc->n;
79  itermax = pastix_data->iparm[IPARM_ITERMAX];
80  eps = pastix_data->dparm[DPARM_EPSILON_REFINEMENT];
81 
82  /* Initialize vectors */
83  gradr = (pastix_complex32_t *)solver.malloc(n * sizeof(pastix_complex32_t));
84  gradp = (pastix_complex32_t *)solver.malloc(n * sizeof(pastix_complex32_t));
85  gradz = (pastix_complex32_t *)solver.malloc(n * sizeof(pastix_complex32_t));
86  grad2 = (pastix_complex32_t *)solver.malloc(n * sizeof(pastix_complex32_t));
87 
88  /* Allocating a vector at half-precision, NULL pointer otherwise */
89  if ( pastix_data->iparm[IPARM_MIXED] )
90  {
91  sgrad = solver.malloc( n * sizeof(pastix_complex32_t) );
92  }
93 
94  clockInit(refine_clk);
95  clockStart(refine_clk);
96 
97  normb = solver.norm( pastix_data, n, b );
98  if ( normb == 0. ) {
99  normb = 1;
100  }
101  normx = solver.norm( pastix_data, n, x );
102 
103  /* Compute r0 = b - A * x */
104  solver.copy( pastix_data, n, b, gradr );
105  if ( normx > 0. ) {
106  solver.spmv( pastix_data, PastixNoTrans, -1., x, 1., gradr );
107  }
108  normr = solver.norm( pastix_data, n, gradr );
109  resid_b = normr / normb;
110 
111  /* z = M^{-1} r */
112  solver.copy( pastix_data, n, gradr, gradz );
113  if ( precond ) {
114  solver.spsv( pastix_data, gradz, sgrad );
115  }
116 
117  /* p = z */
118  solver.copy( pastix_data, n, gradz, gradp );
119 
120  while ((resid_b > eps) && (nb_iter < itermax))
121  {
122  clockStop((refine_clk));
123  t0 = clockGet();
124  nb_iter++;
125 
126  /* grad2 = A * p */
127  solver.spmv( pastix_data, PastixNoTrans, 1.0, gradp, 0., grad2 );
128 
129  /* alpha = <r, z> / <Ap, p> */
130  beta = solver.dot( pastix_data, n, gradr, gradz );
131  alpha = solver.dot( pastix_data, n, grad2, gradp );
132  alpha = beta / alpha;
133 
134  /* x = x + alpha * p */
135  solver.axpy( pastix_data, n, alpha, gradp, x );
136 
137  /* r = r - alpha * A * p */
138  solver.axpy( pastix_data, n, -alpha, grad2, gradr );
139 
140  /* z = M-1 * r */
141  solver.copy( pastix_data, n, gradr, gradz );
142  if ( precond ) {
143  solver.spsv( pastix_data, gradz, sgrad );
144  }
145 
146  /* beta = <r', z> / <r, z> */
147  alpha = solver.dot( pastix_data, n, gradr, gradz );
148  beta = alpha / beta;
149 
150  /* p = z + beta * p */
151  solver.scal( pastix_data, n, beta, gradp );
152  solver.axpy( pastix_data, n, 1., gradz, gradp );
153 
154  normr = solver.norm( pastix_data, n, gradr );
155  resid_b = normr / normb;
156 
157  clockStop((refine_clk));
158  t3 = clockGet();
159  if ( ( pastix_data->iparm[IPARM_VERBOSE] > PastixVerboseNot ) &&
160  ( pastix_data->procnum == 0 ) ) {
161  solver.output_oneiter( t0, t3, resid_b, nb_iter );
162  }
163  t0 = t3;
164  }
165 
166  solver.output_final(pastix_data, resid_b, nb_iter, t3, x, x);
167 
168  solver.free((void*) gradr);
169  solver.free((void*) gradp);
170  solver.free((void*) gradz);
171  solver.free((void*) grad2);
172  solver.free((void*) sgrad);
173 
174  return nb_iter;
175 }
BEGIN_C_DECLS typedef int pastix_int_t
Definition: datatypes.h:51
float _Complex pastix_complex32_t
Definition: datatypes.h:76
double pastix_fixdbl_t
Definition: datatypes.h:65
@ DPARM_EPSILON_REFINEMENT
Definition: api.h:161
@ IPARM_MIXED
Definition: api.h:139
@ IPARM_ITERMAX
Definition: api.h:113
@ IPARM_VERBOSE
Definition: api.h:36
@ PastixNoTrans
Definition: api.h:445
@ PastixVerboseNot
Definition: api.h:220
void c_refine_init(struct c_solver *, pastix_data_t *)
Initiate functions pointers to define basic operations.
pastix_int_t c_grad_smp(pastix_data_t *pastix_data, pastix_rhs_t xp, pastix_rhs_t bp)
Definition: c_refine_grad.c:49
pastix_int_t * iparm
Definition: pastixdata.h:69
double * dparm
Definition: pastixdata.h:70
pastix_bcsc_t * bcsc
Definition: pastixdata.h:101
pastix_int_t steps
Definition: pastixdata.h:72
Main PaStiX data structure.
Definition: pastixdata.h:67
Main PaStiX RHS structure.
Definition: pastixdata.h:150