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