53 struct z_solver solver;
61 pastix_complex64_t *x = (pastix_complex64_t*)(xp->
b);
62 pastix_complex64_t *b = (pastix_complex64_t*)(bp->
b);
63 pastix_complex64_t *gradr;
64 pastix_complex64_t *gradp;
65 pastix_complex64_t *gradz;
66 pastix_complex64_t *grad2;
68 double normb, normx, normr, alpha, beta;
71 memset( &solver, 0,
sizeof(
struct z_solver) );
74 if ( !(pastix_data->
steps & STEP_NUMFACT) ) {
78 n = pastix_data->
bcsc->n;
83 gradr = (pastix_complex64_t *)solver.malloc(n *
sizeof(pastix_complex64_t));
84 gradp = (pastix_complex64_t *)solver.malloc(n *
sizeof(pastix_complex64_t));
85 gradz = (pastix_complex64_t *)solver.malloc(n *
sizeof(pastix_complex64_t));
86 grad2 = (pastix_complex64_t *)solver.malloc(n *
sizeof(pastix_complex64_t));
94 clockInit(refine_clk);
95 clockStart(refine_clk);
97 normb = solver.norm( pastix_data, n, b );
101 normx = solver.norm( pastix_data, n, x );
104 solver.copy( pastix_data, n, b, gradr );
106 solver.spmv( pastix_data,
PastixNoTrans, -1., x, 1., gradr );
108 normr = solver.norm( pastix_data, n, gradr );
109 resid_b = normr / normb;
112 solver.copy( pastix_data, n, gradr, gradz );
114 solver.spsv( pastix_data, gradz, sgrad );
118 solver.copy( pastix_data, n, gradz, gradp );
120 while ((resid_b > eps) && (nb_iter < itermax))
122 clockStop((refine_clk));
127 solver.spmv( pastix_data,
PastixNoTrans, 1.0, gradp, 0., grad2 );
130 beta = solver.dot( pastix_data, n, gradr, gradz );
131 alpha = solver.dot( pastix_data, n, grad2, gradp );
132 alpha = beta / alpha;
135 solver.axpy( pastix_data, n, alpha, gradp, x );
138 solver.axpy( pastix_data, n, -alpha, grad2, gradr );
141 solver.copy( pastix_data, n, gradr, gradz );
143 solver.spsv( pastix_data, gradz, sgrad );
147 alpha = solver.dot( pastix_data, n, gradr, gradz );
151 solver.scal( pastix_data, n, beta, gradp );
152 solver.axpy( pastix_data, n, 1., gradz, gradp );
154 normr = solver.norm( pastix_data, n, gradr );
155 resid_b = normr / normb;
157 clockStop((refine_clk));
160 ( pastix_data->
procnum == 0 ) ) {
161 solver.output_oneiter( t0, t3, resid_b, nb_iter );
166 solver.output_final(pastix_data, resid_b, nb_iter, t3, x, x);
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);
BEGIN_C_DECLS typedef int pastix_int_t
float _Complex pastix_complex32_t
@ DPARM_EPSILON_REFINEMENT
void z_refine_init(struct z_solver *, pastix_data_t *)
Initiate functions pointers to define basic operations.
pastix_int_t z_grad_smp(pastix_data_t *pastix_data, pastix_rhs_t xp, pastix_rhs_t bp)
Main PaStiX data structure.
Main PaStiX RHS structure.