PaStiX Handbook  6.4.0
core_clrdbg.c
Go to the documentation of this file.
1 /**
2  *
3  * @file core_clrdbg.c
4  *
5  * PaStiX low-rank kernel debug routines that may be call within gdb.
6  *
7  * @copyright 2016-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8  * Univ. Bordeaux. All rights reserved.
9  *
10  * @version 6.4.0
11  * @author Gregoire Pichon
12  * @author Mathieu Faverge
13  * @author Pierre Ramet
14  * @date 2024-07-05
15  * @generated from /builds/solverstack/pastix/kernels/core_zlrdbg.c, normal z -> c, Fri Jul 12 15:09:43 2024
16  *
17  **/
18 #include "common.h"
19 #include <cblas.h>
20 #include <lapacke.h>
21 #include "pastix_clrcores.h"
22 #include "c_nan_check.h"
23 
24 /**
25  *******************************************************************************
26  *
27  * @brief Print the svd values of the given matrix.
28  *
29  *******************************************************************************
30  *
31  * @param[in] M
32  * The number of rows of the matrix A.
33  *
34  * @param[in] N
35  * The number of columns of the matrix A.
36  *
37  * @param[in] A
38  * The matrix A to study of size lda-by-N
39  *
40  * @param[in] lda
41  * The leading dimension of the matrix A. lda = max( 1, M )
42  *
43  *******************************************************************************/
44 void
46  pastix_int_t N,
47  const pastix_complex32_t *A,
48  pastix_int_t lda )
49 {
50  pastix_int_t i, ret;
51  pastix_int_t minMN = pastix_imin( M, N );
52  size_t lrwork = 2 * minMN;
53  size_t lzwork = M * N;
55  float *s, *superb;
56 
57  W = malloc( lzwork * sizeof(pastix_complex32_t) + lrwork * sizeof(float) );
58  s = (float*)(W + M*N);
59  superb = s + minMN;
60 
61  ret = LAPACKE_clacpy_work( LAPACK_COL_MAJOR, 'A', M, N, A, lda, W, M );
62  assert( ret == 0 );
63  ret = LAPACKE_cgesvd(LAPACK_COL_MAJOR, 'N', 'N', M, N, W, M, s, NULL, 1, NULL, 1, superb );
64  assert( ret == 0 );
65 
66  for(i=0; i<minMN; i++) {
67  fprintf( stderr, "%e ", s[i] );
68  }
69  fprintf(stderr, "\n");
70 
71  (void)ret;
72  free(W);
73 }
74 
75 /**
76  *******************************************************************************
77  *
78  * @brief Check the orthogonality of the matrix A
79  *
80  *******************************************************************************
81  *
82  * @param[in] M
83  * The number of rows of the matrix A.
84  *
85  * @param[in] N
86  * The number of columns of the matrix A.
87  *
88  * @param[in] A
89  * The matrix A to study of size lda-by-N
90  *
91  * @param[in] lda
92  * The leading dimension of the matrix A. lda = max( 1, M )
93  *
94  *******************************************************************************
95  *
96  * @retval 0 if the matrix A is orthogonal
97  * @retval 1 if the matrix A is not orthogonal
98  *
99  *******************************************************************************/
100 int
102  pastix_int_t N,
103  const pastix_complex32_t *A,
104  pastix_int_t lda )
105 {
106  pastix_complex32_t *Id;
107  float alpha, beta;
108  float normQ, res;
109  pastix_int_t info_ortho, ret;
110  pastix_int_t minMN = pastix_imin(M, N);
111  pastix_int_t maxMN = pastix_imax(M, N);
112  float eps = LAPACKE_slamch_work('e');
113 
114  alpha = 1.0;
115  beta = -1.0;
116 
117  /* Build the identity matrix */
118  Id = malloc( minMN * minMN * sizeof(pastix_complex32_t) );
119  ret = LAPACKE_claset_work( LAPACK_COL_MAJOR, 'A', minMN, minMN,
120  0., 1., Id, minMN );
121  assert( ret == 0 );
122 
123  if (M > N) {
124  /* Perform Id - Q'Q */
125  cblas_cherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, alpha, A, lda, beta, Id, minMN);
126  }
127  else {
128  /* Perform Id - QQ' */
129  cblas_cherk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, A, lda, beta, Id, minMN);
130  }
131 
132  normQ = LAPACKE_clanhe_work( LAPACK_COL_MAJOR, 'M', 'U', minMN, Id, minMN, NULL );
133  res = normQ / (maxMN * eps);
134 
135  if ( isnan(res) || isinf(res) || (res > 60.0) ) {
136  fprintf(stderr, "Check Orthogonality: || I - Q*Q' || = %e, ||Id-Q'*Q||_oo / (N*eps) = %e : \n",
137  normQ, res );
138  info_ortho = 1;
139  }
140  else {
141  info_ortho = 0;
142  }
143 
144  free(Id);
145  (void)ret;
146  return info_ortho;
147 }
148 
149 /**
150  *******************************************************************************
151  *
152  * @brief Check the orthogonality of the matrix A relatively to the matrix B
153  *
154  * Check that A^t B = 0
155  *
156  *******************************************************************************
157  *
158  * @param[in] M
159  * The number of rows of the matrix A.
160  *
161  * @param[in] NA
162  * The number of columns of the matrix A.
163  *
164  * @param[in] NB
165  * The number of columns of the matrix B.
166  *
167  * @param[in] A
168  * The matrix A to study of size lda-by-NA
169  *
170  * @param[in] lda
171  * The leading dimension of the matrix A. lda = max( 1, M )
172  *
173  * @param[in] B
174  * The matrix B to study of size ldb-by-NB
175  *
176  * @param[in] ldb
177  * The leading dimension of the matrix B. ldb = max( 1, M )
178  *
179  *******************************************************************************
180  *
181  * @retval 0 if the matrices A and B are orthogonal
182  * @retval 1 if the matrices A anb B are not orthogonal
183  *
184  *******************************************************************************/
185 int
187  pastix_int_t NA,
188  pastix_int_t NB,
189  const pastix_complex32_t *A,
190  pastix_int_t lda,
191  const pastix_complex32_t *B,
192  pastix_int_t ldb )
193 {
194  pastix_complex32_t *Zero;
195  float norm, res;
196  pastix_int_t info_ortho, ret;
197  float eps = LAPACKE_slamch_work('e');
198  pastix_complex32_t cone = 1.0;
199  pastix_complex32_t czero = 0.0;
200 
201  /* Build the null matrix */
202  Zero = malloc( NA * NB * sizeof(pastix_complex32_t) );
203  ret = LAPACKE_claset_work( LAPACK_COL_MAJOR, 'A', NA, NB,
204  0., 0., Zero, NA );
205  assert( ret == 0 );
206 
207  cblas_cgemm(CblasColMajor, CblasConjTrans, CblasNoTrans,
208  NA, NB, M,
209  CBLAS_SADDR(cone), A, lda,
210  B, ldb,
211  CBLAS_SADDR(czero), Zero, NA);
212 
213  norm = LAPACKE_clange_work( LAPACK_COL_MAJOR, 'M', NA, NB, Zero, NA, NULL );
214  res = norm / (M * eps);
215 
216  if ( isnan(res) || isinf(res) || (res > 60.0) ) {
217  fprintf(stderr, "Check Orthogonality: || A' B || = %e, || A' B ||_oo / (M*eps) = %e : \n",
218  norm, res );
219  info_ortho = 1;
220  }
221  else {
222  info_ortho = 0;
223  }
224 
225  free(Zero);
226  (void)ret;
227  return info_ortho;
228 }
Manage nancheck for lowrank kernels. This header describes all the LAPACKE functions used for low-ran...
BEGIN_C_DECLS typedef int pastix_int_t
Definition: datatypes.h:51
float _Complex pastix_complex32_t
Definition: datatypes.h:76
int core_clrdbg_check_orthogonality_AB(pastix_int_t M, pastix_int_t NA, pastix_int_t NB, const pastix_complex32_t *A, pastix_int_t lda, const pastix_complex32_t *B, pastix_int_t ldb)
Check the orthogonality of the matrix A relatively to the matrix B.
Definition: core_clrdbg.c:186
int core_clrdbg_check_orthogonality(pastix_int_t M, pastix_int_t N, const pastix_complex32_t *A, pastix_int_t lda)
Check the orthogonality of the matrix A.
Definition: core_clrdbg.c:101
void core_clrdbg_printsvd(pastix_int_t M, pastix_int_t N, const pastix_complex32_t *A, pastix_int_t lda)
Print the svd values of the given matrix.
Definition: core_clrdbg.c:45