SpM Handbook 1.2.4
Loading...
Searching...
No Matches
c_spm_genrhs.c
Go to the documentation of this file.
1/**
2 *
3 * @file c_spm_genrhs.c
4 *
5 * SParse Matrix package right hand side generators.
6 *
7 * @copyright 2016-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8 * Univ. Bordeaux. All rights reserved.
9 *
10 * @version 1.2.4
11 * @author Mathieu Faverge
12 * @author Tony Delarue
13 * @date 2024-06-08
14 *
15 * @generated from /builds/2mk6rsew/0/fpruvost/spm/src/z_spm_genrhs.c, normal z -> c, Fri Nov 29 11:34:28 2024
16 **/
17#include "common.h"
18#include <cblas.h>
19#include <lapacke.h>
20
21#ifndef DOXYGEN_SHOULD_SKIP_THIS
22static spm_complex32_t mcone = (spm_complex32_t)-1.;
23#endif /* DOXYGEN_SHOULD_SKIP_THIS */
24
25/**
26 *******************************************************************************
27 *
28 * @ingroup spm_dev_rhs
29 *
30 * @brief Generate nrhs right hand side vectors associated to a given
31 * matrix to test a problem with a solver.
32 *
33 *******************************************************************************
34 *
35 * @param[in] type
36 * Defines how to compute the vector b.
37 * - SpmRhsOne: b is computed such that x = 1 [ + I ]
38 * - SpmRhsI: b is computed such that x = i [ + i * I ]
39 * - SpmRhsRndX: b is computed by matrix-vector product, such that
40 * is a random vector in the range [-0.5, 0.5]
41 * - SpmRhsRndB: b is computed randomly and x is not computed.
42 *
43 * @param[in] nrhs
44 * Defines the number of right hand side that must be generated.
45 *
46 * @param[in] spm
47 * The sparse matrix uses to generate the right hand side, and the
48 * solution of the full problem.
49 *
50 * @param[out] x
51 * On exit, if x != NULL, then the x vector(s) generated to compute b
52 * is returned. Must be of size at least ldx * spm->n.
53 *
54 * @param[in] ldx
55 * Defines the leading dimension of x when multiple right hand sides
56 * are available. ldx >= spm->n.
57 *
58 * @param[inout] b
59 * b must be an allocated matrix of size at least ldb * nrhs.
60 * On exit, b is initialized as defined by the type parameter.
61 *
62 * @param[in] ldb
63 * Defines the leading dimension of b when multiple right hand sides
64 * are available. ldb >= spm->n.
65 *
66 *******************************************************************************
67 *
68 * @retval SPM_SUCCESS if the b vector has been computed succesfully,
69 * @retval SPM_ERR_BADPARAMETER otherwise.
70 *
71 *******************************************************************************/
72int
74 int nrhs,
75 const spmatrix_t *spm,
76 void *x,
77 int ldx,
78 void *b,
79 int ldb )
80{
84 int rc;
85
86 if( ( spm == NULL ) ||
87 ( spm->values == NULL ) ) {
89 }
90
91 if( spm->gN <= 0 ) {
93 }
94
95 if( nrhs <= 0 ) {
97 }
98
99 if( (xptr != NULL) && (nrhs > 1) && (ldx < spm->nexp) ) {
101 }
102
103 if( b == NULL ) {
105 }
106
107 if( ( nrhs > 1 ) && ( ldb < spm->nexp ) ) {
109 }
110
111 if ( nrhs == 1 ) {
112 ldb = spm->nexp;
113 ldx = spm->nexp;
114 }
115
116 /* If random b, we do it and exit */
117 if ( type == SpmRhsRndB ) {
118 /* Compute the spm norm to scale the b vector */
120 if ( norm == 0. ) {
121 norm = 1.;
122 }
123 c_spmGenMat( type, nrhs, spm, &norm, 24356, bptr, ldb );
124
125 return SPM_SUCCESS;
126 }
127
128 if ( (type == SpmRhsOne ) ||
129 (type == SpmRhsI ) ||
130 (type == SpmRhsRndX ) )
131 {
132 if ( xptr == NULL ) {
133 ldx = spm->nexp;
134 xptr = malloc( ldx * nrhs * sizeof(spm_complex32_t) );
135 }
136
137 c_spmGenMat( type, nrhs, spm, &alpha, 24356, xptr, ldx );
138
139 /* Compute B */
140 rc = spm_cspmm( SpmLeft, SpmNoTrans, SpmNoTrans, nrhs, alpha, spm, xptr, ldx, 0., bptr, ldb );
141
142 if ( x == NULL ) {
143 free(xptr);
144 }
145 return rc;
146 }
147
148 fprintf(stderr, "c_spmGenRHS: Generator not implemented yet\n");
149
150 return SPM_SUCCESS;
151}
152
153/**
154 *******************************************************************************
155 *
156 * @ingroup spm_dev_rhs
157 *
158 * @brief Check the backward error, and the forward error if x0 is provided.
159 *
160 *******************************************************************************
161 *
162 * @param[in] eps
163 * The epsilon threshold used for the refinement step. -1. to use the
164 * machine precision.
165 *
166 * @param[in] nrhs
167 * Defines the number of right hand side that must be generated.
168 *
169 * @param[in] spm
170 * The sparse matrix uses to generate the right hand side, and the
171 * solution of the full problem.
172 *
173 * @param[inout] x0
174 * If x0 != NULL, the forward error is computed.
175 * On exit, x0 stores (x0-x)
176 *
177 * @param[in] ldx0
178 * Defines the leading dimension of x0 when multiple right hand sides
179 * are available. ldx0 >= spm->n.
180 *
181 * @param[inout] b
182 * b is a matrix of size at least ldb * nrhs.
183 * On exit, b stores Ax-b.
184 *
185 * @param[in] ldb
186 * Defines the leading dimension of b when multiple right hand sides
187 * are available. ldb >= spm->n.
188 *
189 * @param[in] x
190 * Contains the solution computed by the solver.
191 *
192 * @param[in] ldx
193 * Defines the leading dimension of x when multiple right hand sides
194 * are available. ldx >= spm->n.
195 *
196 *******************************************************************************
197 *
198 * @retval SPM_SUCCESS if the tests are succesfull
199 * @retval 1, if one of the test failed
200 *
201 *******************************************************************************/
202int
204 const spmatrix_t *spm,
205 void *x0, int ldx0,
206 void *b, int ldb,
207 const void *x, int ldx )
208{
209 const spm_complex32_t *zx = (const spm_complex32_t *)x;
210 spm_complex32_t *zx0 = (spm_complex32_t *)x0;
212 float *nb2 = malloc( nrhs * sizeof(float) );
213 float normA, normB, normX, normR, normR2;
214 float backward, forward;
215 int failure = 0;
216 int i;
217
218 if ( eps == -1. ) {
219 eps = LAPACKE_slamch('e');
220 }
221
222 /**
223 * Compute the starting norms
224 */
225 normA = spmNorm( SpmOneNorm, spm );
226
227 normB = 0.;
228 normX = 0.;
229 for( i=0; i<nrhs; i++ ) {
230 float norm;
231
232 norm = c_spmNormMat( SpmInfNorm, spm, 1, zb + i * ldb, ldb );
233 normB = ( norm > normB ) ? norm : normB;
234 norm = c_spmNormMat( SpmInfNorm, spm, 1, zx + i * ldx, ldx );
235 normX = ( norm > normX ) ? norm : normX;
236
237 nb2[i] = c_spmNormMat( SpmFrobeniusNorm, spm, 1, zb + i * ldb, ldb );
238 if ( nb2[i] == 0. ) { nb2[i] = 1.; }
239 }
240 fprintf( stdout,
241 " || A ||_1 %e\n"
242 " max(|| b_i ||_oo) %e\n"
243 " max(|| x_i ||_oo) %e\n",
244 normA, normB, normX );
245
246 /**
247 * Compute r = b - A * x
248 */
249 spm_cspmm( SpmLeft, SpmNoTrans, SpmNoTrans, nrhs, -1., spm, x, ldx, 1., b, ldb );
250
251 normR = 0.;
252 normR2 = 0.;
253 backward = 0.;
254 failure = 0;
255
256 for( i=0; i<nrhs; i++ ) {
257 float nx = c_spmNormMat( SpmOneNorm, spm, 1, zx + i * ldx, ldx );
258 float nr = c_spmNormMat( SpmOneNorm, spm, 1, zb + i * ldb, ldb );
259 float nr2 = c_spmNormMat( SpmFrobeniusNorm, spm, 1, zb + i * ldb, ldb ) / nb2[i];
260 float back = ( nr / eps );
261 int fail = 0;
262
263 if ( normA > 0. ) {
264 nr = nr / normA;
265 }
266 if ( nx > 0. ) {
267 nr = nr / nx;
268 }
269
270 normR = (nr > normR ) ? nr : normR;
271 normR2 = (nr2 > normR2 ) ? nr2 : normR2;
272 backward = (back > backward) ? back : backward;
273
274 fail = isnan(nr) || isinf(nr) || isnan(back) || isinf(back) || (back > 1.e2);
275 if ( fail ) {
276 fprintf( stdout,
277 " || b_%d - A x_%d ||_2 / || b_%d ||_2 %e\n"
278 " || b_%d - A x_%d ||_1 %e\n"
279 " || b_%d - A x_%d ||_1 / (||A||_1 * ||x_%d||_oo * eps) %e (FAILED)\n",
280 i, i, i, nr2,
281 i, i, nr,
282 i, i, i, back );
283 }
284
285 failure = failure || fail;
286 }
287
288 fprintf( stdout,
289 " max(|| b_i - A x_i ||_2 / || b_i ||_2) %e\n"
290 " max(|| b_i - A x_i ||_1) %e\n"
291 " max(|| b_i - A x_i ||_1 / (||A||_1 * ||x_i||_oo * eps)) %e (%s)\n",
292 normR2, normR, backward,
293 failure ? "FAILED" : "SUCCESS" );
294
295 free(nb2);
296
297 /**
298 * Compute r = x0 - x
299 */
300 if ( x0 != NULL ) {
301 float normX0;
302 float forw, nr, nx, nx0;
303 int fail;
304
305 forward = 0.;
306 normR = 0.;
307 normX0 = 0.;
308 failure = 0;
309
310 for( i=0; i<nrhs; i++, zx += ldx, zx0 += ldx0 ) {
311
312 nx0 = c_spmNormMat( SpmInfNorm, spm, 1, zx0, ldx0 );
313 nx = c_spmNormMat( SpmInfNorm, spm, 1, zx, ldx );
314
315 cblas_caxpy( spm->nexp, CBLAS_SADDR(mcone),
316 zx, 1, zx0, 1);
317
318 nr = c_spmNormMat( SpmInfNorm, spm, 1, zx0, ldx0 );
319
320 forw = nr / eps;
321 if ( nx0 > 0. ) { forw = forw / nx0; }
322
323 normX0 = ( nx > normX0 ) ? nx : normX0;
324 normR = ( nr > normR ) ? nr : normR;
325 forward = ( forw > forward ) ? forw : forward;
326
327 fail = isnan(nx) || isinf(nx) || isnan(forw) || isinf(forw) || (forw > 1.e2);
328 if ( fail ) {
329 fprintf( stdout,
330 " || x_%d ||_oo %e\n"
331 " || x0_%d - x_%d ||_oo %e\n"
332 " || x0_%d - x_%d ||_oo / (||x0_%d||_oo * eps) %e (FAILED)\n",
333 i, nx,
334 i, i, nr,
335 i, i, i, forw );
336 }
337
338 failure = failure || fail;
339 }
340
341 fprintf( stdout,
342 " max(|| x_i ||_oo) %e\n"
343 " max(|| x0_i - x_i ||_oo) %e\n"
344 " max(|| x0_i - x_i ||_oo / || x0_i ||_oo) %e (%s)\n",
345 normX0, normR, forward,
346 failure ? "FAILED" : "SUCCESS" );
347 }
348
349 fflush( stdout );
350
351 return - failure;
352}
double spm_fixdbl_t
Double datatype that is not converted through precision generator functions.
Definition datatypes.h:150
float _Complex spm_complex32_t
Definition datatypes.h:161
#define CBLAS_SADDR(_a_)
Macro to get the address if need in the cblas calls.
Definition const.h:36
enum spm_rhstype_e spm_rhstype_t
How to generate RHS.
@ SPM_ERR_BADPARAMETER
Definition const.h:89
@ SPM_SUCCESS
Definition const.h:82
@ SpmNoTrans
Definition const.h:155
@ SpmFrobeniusNorm
Definition const.h:200
@ SpmInfNorm
Definition const.h:201
@ SpmOneNorm
Definition const.h:199
@ SpmLeft
Definition const.h:191
int spm_cspmm(spm_side_t side, spm_trans_t transA, spm_trans_t transB, spm_int_t K, spm_complex32_t alpha, const spmatrix_t *A, const spm_complex32_t *B, spm_int_t ldb, spm_complex32_t beta, spm_complex32_t *C, spm_int_t ldc)
Compute a matrix-matrix product.
float c_spmNormMat(spm_normtype_t ntype, const spmatrix_t *spm, spm_int_t n, const spm_complex32_t *A, spm_int_t lda)
Compute the norm of a dense matrix that follows the distribution of an spm matrix.
float c_spmNorm(spm_normtype_t ntype, const spmatrix_t *spm)
Compute the norm of an spm matrix.
int c_spmGenRHS(spm_rhstype_t type, int nrhs, const spmatrix_t *spm, void *x, int ldx, void *b, int ldb)
Generate nrhs right hand side vectors associated to a given matrix to test a problem with a solver.
int c_spmGenMat(spm_rhstype_t type, int nrhs, const spmatrix_t *spm, void *alpha, unsigned long long int seed, void *A, int lda)
Generate nrhs right hand side vectors associated to a given matrix to test a problem with a solver.
int c_spmCheckAxb(spm_fixdbl_t eps, int nrhs, const spmatrix_t *spm, void *x0, int ldx0, void *b, int ldb, const void *x, int ldx)
Check the backward error, and the forward error if x0 is provided.
spm_int_t nexp
Definition spm.h:78
void * values
Definition spm.h:92
spm_int_t gN
Definition spm.h:72
double spmNorm(spm_normtype_t ntype, const spmatrix_t *spm)
Compute the norm of the spm.
Definition spm.c:514
The sparse matrix data structure.
Definition spm.h:64