PaStiX Handbook  6.3.1
s_refine_pivot.c
Go to the documentation of this file.
1 /**
2  *
3  * @file s_refine_pivot.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.1
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 2023-07-21
18  * @generated from /builds/solverstack/pastix/refinement/z_refine_pivot.c, normal z -> s, Thu Nov 23 09:56:16 2023
19  *
20  **/
21 #include "common.h"
22 #include "s_refine_functions.h"
23 
24 /**
25  *******************************************************************************
26  *
27  * @ingroup pastix_refine
28  *
29  * s_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 s_pivot_smp( pastix_data_t *pastix_data,
49  pastix_rhs_t xp,
50  pastix_rhs_t bp )
51 {
52  struct s_solver solver;
53  Clock t0, t3, refine_clk;
54  pastix_int_t n, itermax;
55  float *x = (float*)(xp->b);
56  float *b = (float*)(bp->b);
57  float *r, *dx;
58  float *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 s_solver) );
65  s_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 = (float *)solver.malloc(n * sizeof(float));
81  dx = (float *)solver.malloc(n * sizeof(float));
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(float) );
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
@ 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 s_refine_init(struct s_solver *, pastix_data_t *)
Initiate functions pointers to define basic operations.
pastix_int_t s_pivot_smp(pastix_data_t *pastix_data, pastix_rhs_t xp, pastix_rhs_t bp)
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