1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
/*
-- MAGMA (version 2.9.0) --
Univ. of Tennessee, Knoxville
Univ. of California, Berkeley
Univ. of Colorado, Denver
@date January 2025
@author Mark Gates
@generated from testing/magma_zgesvd_check.cpp, normal z -> d, Wed Jan 22 14:40:18 2025
*/
#include "magma_v2.h"
#include "magma_lapack.h"
#include "testings.h"
#define REAL
/**
Check the results following the LAPACK's [zcds]drvbd routine.
A is factored as A = U diag(S) VT and the following 4 tests computed:
(1) | A - U diag(S) VT | / ( |A| max(m,n) )
(2) | I - U^H U | / ( m )
(3) | I - VT VT^H | / ( n )
(4) S contains min(m,n) nonnegative values in decreasing order.
(Return 0 if true, 1 if false.)
If check is false, skips (1) - (3), but always does (4).
********************************************************************/
extern "C"
void check_dgesvd(
magma_int_t check,
magma_vec_t jobu,
magma_vec_t jobv,
magma_int_t m, magma_int_t n,
double *A, magma_int_t lda,
double *S,
double *U, magma_int_t ldu,
double *VT, magma_int_t ldv,
double result[4] )
{
double unused[1];
const magma_int_t izero = 0;
double eps = lapackf77_dlamch( "E" );
if ( jobu == MagmaNoVec ) {
U = NULL;
}
if ( jobv == MagmaNoVec ) {
VT = NULL;
}
// -1 indicates check not done
result[0] = -1;
result[1] = -1;
result[2] = -1;
result[3] = -1;
magma_int_t min_mn = min(m, n);
magma_int_t n_u = (jobu == MagmaAllVec ? m : min_mn);
magma_int_t m_vt = (jobv == MagmaAllVec ? n : min_mn);
assert( lda >= m );
assert( ldu >= m );
assert( ldv >= m_vt );
if ( check ) {
// dbdt01 needs m+n
// dort01 prefers n*(n+1) to check U; m*(m+1) to check V
magma_int_t lwork_err = m+n;
if ( U != NULL ) {
lwork_err = max( lwork_err, n_u*(n_u+1) );
}
if ( VT != NULL ) {
lwork_err = max( lwork_err, m_vt*(m_vt+1) );
}
double *work_err;
TESTING_CHECK( magma_dmalloc_cpu( &work_err, lwork_err ));
// dbdt01 and dort01 need max(m,n), depending
double *rwork_err;
TESTING_CHECK( magma_dmalloc_cpu( &rwork_err, max(m,n) ));
if ( U != NULL && VT != NULL ) {
// since KD=0 (3rd arg), E is not referenced so pass unused (9th arg)
lapackf77_dbdt01( &m, &n, &izero, A, &lda,
U, &ldu, S, unused, VT, &ldv,
work_err,
#ifdef COMPLEX
rwork_err,
#endif
&result[0] );
}
if ( U != NULL ) {
lapackf77_dort01( "Columns", &m, &n_u, U, &ldu, work_err, &lwork_err,
#ifdef COMPLEX
rwork_err,
#endif
&result[1] );
}
if ( VT != NULL ) {
lapackf77_dort01( "Rows", &m_vt, &n, VT, &ldv, work_err, &lwork_err,
#ifdef COMPLEX
rwork_err,
#endif
&result[2] );
}
result[0] *= eps;
result[1] *= eps;
result[2] *= eps;
magma_free_cpu( work_err );
magma_free_cpu( rwork_err );
}
// check S is sorted
result[3] = 0.;
for (int j=0; j < min_mn-1; j++) {
if ( S[j] < S[j+1] )
result[3] = 1.;
if ( S[j] < 0. )
result[3] = 1.;
}
if (min_mn > 1 && S[min_mn-1] < 0.) {
result[3] = 1.;
}
}
|
