PaStiX Handbook  6.4.0
c_refine_pivot.c
Go to the documentation of this file.
1 /**
2  *
3  * @file c_refine_pivot.c
4  *
5  * PaStiX refinement functions implementations.
6  *
7  * @copyright 2015-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8  * Univ. Bordeaux. All rights reserved.
9  *
10  * @version 6.4.0
11  * @author Mathieu Faverge
12  * @author Pierre Ramet
13  * @author Xavier Lacoste
14  * @author Gregoire Pichon
15  * @author Theophile Terraz
16  * @author Vincent Bridonneau
17  * @date 2024-07-05
18  * @generated from /builds/solverstack/pastix/refinement/z_refine_pivot.c, normal z -> c, Fri Jul 12 15:10:19 2024
19  *
20  **/
21 #include "common.h"
22 #include "c_refine_functions.h"
23 
24 /**
25  *******************************************************************************
26  *
27  * @ingroup pastix_refine
28  *
29  * c_pivot_smp - Refine the solution using static pivoting method.
30  *
31  *******************************************************************************
32  *
33  * @param[in] pastix_data
34  * The PaStiX data structure that describes the solver instance.
35  *
36  * @param[out] xp
37  * The solution vector.
38  *
39  * @param[in] bp
40  * The right hand side member (only one).
41  *
42  *******************************************************************************
43  *
44  * @return Number of iterations
45  *
46  *******************************************************************************/
48 c_pivot_smp( pastix_data_t *pastix_data,
49  pastix_rhs_t xp,
50  pastix_rhs_t bp )
51 {
52  struct c_solver solver;
53  Clock t0, t3, refine_clk;
54  pastix_int_t n, itermax;
57  pastix_complex32_t *r, *dx;
58  pastix_complex32_t *sb = NULL;
59  float eps, normb, normr;
60  float berr, last_berr;
61  pastix_int_t iter = 0;
62  int flag = 1;
63 
64  memset( &solver, 0, sizeof(struct c_solver) );
65  c_refine_init( &solver, pastix_data );
66 
67  if ( !(pastix_data->steps & STEP_NUMFACT) ) {
68  fprintf(stderr, "pastix_task_refine: Simple refinement cannot be applied without preconditionner\n" );
69  return -1;
70  }
71 
72  n = pastix_data->bcsc->n;
73  itermax = pastix_data->iparm[IPARM_ITERMAX];
74  eps = pastix_data->dparm[DPARM_EPSILON_REFINEMENT];
75 
76  if (pastix_data->iparm[IPARM_VERBOSE] > PastixVerboseNot)
77  {
78  fprintf(stdout, OUT_ITERREFINE_PIVOT);
79  }
80  r = (pastix_complex32_t *)solver.malloc(n * sizeof(pastix_complex32_t));
81  dx = (pastix_complex32_t *)solver.malloc(n * sizeof(pastix_complex32_t));
82 
83  clockInit(refine_clk);
84  clockStart(refine_clk);
85 
86  normb = solver.norm( pastix_data, n, b );
87  if ( normb == 0. ) {
88  normb = 1;
89  }
90 
91  /* Allocating a vector at half-precision, NULL pointer otherwise */
92  if ( pastix_data->iparm[IPARM_MIXED] )
93  {
94  sb = solver.malloc( n * sizeof(pastix_complex32_t) );
95  }
96 
97  t0 = clockGet();
98  while(flag)
99  {
100  /* Compute r = b - A * x */
101  solver.copy( pastix_data, n, b, r );
102  solver.spmv( pastix_data, PastixNoTrans, -1., x, 1., r );
103 
104  /*
105  * berr should be equal to the componentwise backward error in the literature:
106  * max( r_i / ( |A| |x| + |b| )_i )
107  * For simplicity, we replace it by ||r||_f / ||b||_f which may not be
108  * as good as the previous one.
109  */
110  normr = solver.norm( pastix_data, n, r );
111  berr = normr / normb;
112 
113  /* Force the first error */
114  if ( iter == 0 ) {
115  last_berr = 3 * berr;
116  }
117  else {
118  t3 = clockGet();
119  if ( ( pastix_data->iparm[IPARM_VERBOSE] > PastixVerboseNot ) &&
120  ( pastix_data->procnum == 0 ) ) {
121  solver.output_oneiter( t0, t3, berr, iter );
122  }
123  t0 = clockGet();
124  }
125 
126  if ( (iter < itermax) &&
127  (berr > eps) &&
128  (berr <= (last_berr / 2.)) )
129  {
130  t3 = clockGet();
131 
132  /* Solve A dx = r */
133  solver.copy( pastix_data, n, r, dx );
134  solver.spsv( pastix_data, dx, sb );
135 
136  /* Accumulate the solution: x = x + dx */
137  solver.axpy( pastix_data, n, 1.0, dx, x );
138 
139  last_berr = berr;
140  }
141  else
142  {
143  flag = 0;
144  }
145 
146  iter++;
147  }
148  clockStop(refine_clk);
149 
150  solver.output_final( pastix_data, berr, iter, refine_clk, x, x );
151 
152  solver.free( r );
153  solver.free( dx );
154  solver.free( sb );
155 
156  return iter;
157 }
BEGIN_C_DECLS typedef int pastix_int_t
Definition: datatypes.h:51
float _Complex pastix_complex32_t
Definition: datatypes.h:76
@ 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_pivot_smp(pastix_data_t *pastix_data, pastix_rhs_t xp, pastix_rhs_t bp)
pastix_int_t * iparm
Definition: pastixdata.h:70
double * dparm
Definition: pastixdata.h:71
pastix_bcsc_t * bcsc
Definition: pastixdata.h:102
pastix_int_t steps
Definition: pastixdata.h:73
Main PaStiX data structure.
Definition: pastixdata.h:68
Main PaStiX RHS structure.
Definition: pastixdata.h:155