1// Copyright ©2015 The Gonum Authors. All rights reserved. 2// Use of this source code is governed by a BSD-style 3// license that can be found in the LICENSE file. 4 5package testlapack 6 7import ( 8 "fmt" 9 "math" 10 "sort" 11 "testing" 12 13 "golang.org/x/exp/rand" 14 15 "gonum.org/v1/gonum/blas" 16 "gonum.org/v1/gonum/blas/blas64" 17 "gonum.org/v1/gonum/floats" 18 "gonum.org/v1/gonum/lapack" 19) 20 21type Dgesvder interface { 22 Dgesvd(jobU, jobVT lapack.SVDJob, m, n int, a []float64, lda int, s, u []float64, ldu int, vt []float64, ldvt int, work []float64, lwork int) (ok bool) 23} 24 25func DgesvdTest(t *testing.T, impl Dgesvder, tol float64) { 26 for _, m := range []int{0, 1, 2, 3, 4, 5, 10, 150, 300} { 27 for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 150} { 28 for _, mtype := range []int{1, 2, 3, 4, 5} { 29 dgesvdTest(t, impl, m, n, mtype, tol) 30 } 31 } 32 } 33} 34 35// dgesvdTest tests a Dgesvd implementation on an m×n matrix A generated 36// according to mtype as: 37// - the zero matrix if mtype == 1, 38// - the identity matrix if mtype == 2, 39// - a random matrix with a given condition number and singular values if mtype == 3, 4, or 5. 40// It first computes the full SVD A = U*Sigma*V^T and checks that 41// - U has orthonormal columns, and V^T has orthonormal rows, 42// - U*Sigma*V^T multiply back to A, 43// - the singular values are non-negative and sorted in decreasing order. 44// Then all combinations of partial SVD results are computed and checked whether 45// they match the full SVD result. 46func dgesvdTest(t *testing.T, impl Dgesvder, m, n, mtype int, tol float64) { 47 rnd := rand.New(rand.NewSource(1)) 48 49 // Use a fixed leading dimension to reduce testing time. 50 lda := n + 3 51 ldu := m + 5 52 ldvt := n + 7 53 54 minmn := min(m, n) 55 56 // Allocate A and fill it with random values. The in-range elements will 57 // be overwritten below according to mtype. 58 a := make([]float64, m*lda) 59 for i := range a { 60 a[i] = rnd.NormFloat64() 61 } 62 63 var aNorm float64 64 switch mtype { 65 default: 66 panic("unknown test matrix type") 67 case 1: 68 // Zero matrix. 69 for i := 0; i < m; i++ { 70 for j := 0; j < n; j++ { 71 a[i*lda+j] = 0 72 } 73 } 74 aNorm = 0 75 case 2: 76 // Identity matrix. 77 for i := 0; i < m; i++ { 78 for j := 0; j < n; j++ { 79 if i == j { 80 a[i*lda+i] = 1 81 } else { 82 a[i*lda+j] = 0 83 } 84 } 85 } 86 aNorm = 1 87 case 3, 4, 5: 88 // Scaled random matrix. 89 // Generate singular values. 90 s := make([]float64, minmn) 91 Dlatm1(s, 92 4, // s[i] = 1 - i*(1-1/cond)/(minmn-1) 93 float64(max(1, minmn)), // where cond = max(1,minmn) 94 false, // signs of s[i] are not randomly flipped 95 1, rnd) // random numbers are drawn uniformly from [0,1) 96 // Decide scale factor for the singular values based on the matrix type. 97 ulp := dlamchP 98 unfl := dlamchS 99 ovfl := 1 / unfl 100 aNorm = 1 101 if mtype == 4 { 102 aNorm = unfl / ulp 103 } 104 if mtype == 5 { 105 aNorm = ovfl * ulp 106 } 107 // Scale singular values so that the maximum singular value is 108 // equal to aNorm (we know that the singular values are 109 // generated above to be spread linearly between 1/cond and 1). 110 floats.Scale(aNorm, s) 111 // Generate A by multiplying S by random orthogonal matrices 112 // from left and right. 113 Dlagge(m, n, max(0, m-1), max(0, n-1), s, a, lda, rnd, make([]float64, m+n)) 114 } 115 aCopy := make([]float64, len(a)) 116 copy(aCopy, a) 117 118 for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} { 119 // Restore A because Dgesvd overwrites it. 120 copy(a, aCopy) 121 122 // Allocate slices that will be used below to store the results of full 123 // SVD and fill them. 124 uAll := make([]float64, m*ldu) 125 for i := range uAll { 126 uAll[i] = rnd.NormFloat64() 127 } 128 vtAll := make([]float64, n*ldvt) 129 for i := range vtAll { 130 vtAll[i] = rnd.NormFloat64() 131 } 132 sAll := make([]float64, min(m, n)) 133 for i := range sAll { 134 sAll[i] = math.NaN() 135 } 136 137 prefix := fmt.Sprintf("m=%v,n=%v,work=%v,mtype=%v", m, n, wl, mtype) 138 139 // Determine workspace size based on wl. 140 minwork := max(1, max(5*min(m, n), 3*min(m, n)+max(m, n))) 141 var lwork int 142 switch wl { 143 case minimumWork: 144 lwork = minwork 145 case mediumWork: 146 work := make([]float64, 1) 147 impl.Dgesvd(lapack.SVDAll, lapack.SVDAll, m, n, a, lda, sAll, uAll, ldu, vtAll, ldvt, work, -1) 148 lwork = (int(work[0]) + minwork) / 2 149 case optimumWork: 150 work := make([]float64, 1) 151 impl.Dgesvd(lapack.SVDAll, lapack.SVDAll, m, n, a, lda, sAll, uAll, ldu, vtAll, ldvt, work, -1) 152 lwork = int(work[0]) 153 } 154 work := make([]float64, max(1, lwork)) 155 for i := range work { 156 work[i] = math.NaN() 157 } 158 159 // Compute the full SVD which will be used later for checking the partial results. 160 ok := impl.Dgesvd(lapack.SVDAll, lapack.SVDAll, m, n, a, lda, sAll, uAll, ldu, vtAll, ldvt, work, len(work)) 161 if !ok { 162 t.Fatalf("Case %v: unexpected failure in full SVD", prefix) 163 } 164 165 // Check that uAll, sAll, and vtAll multiply back to A by computing a residual 166 // |A - U*S*VT| / (n*aNorm) 167 if resid := svdFullResidual(m, n, aNorm, aCopy, lda, uAll, ldu, sAll, vtAll, ldvt); resid > tol { 168 t.Errorf("Case %v: original matrix not recovered for full SVD, |A - U*D*VT|=%v", prefix, resid) 169 } 170 if minmn > 0 { 171 // Check that uAll is orthogonal. 172 if !hasOrthonormalColumns(blas64.General{Rows: m, Cols: m, Data: uAll, Stride: ldu}) { 173 t.Errorf("Case %v: UAll is not orthogonal", prefix) 174 } 175 // Check that vtAll is orthogonal. 176 if !hasOrthonormalRows(blas64.General{Rows: n, Cols: n, Data: vtAll, Stride: ldvt}) { 177 t.Errorf("Case %v: VTAll is not orthogonal", prefix) 178 } 179 } 180 // Check that singular values are decreasing. 181 if !sort.IsSorted(sort.Reverse(sort.Float64Slice(sAll))) { 182 t.Errorf("Case %v: singular values from full SVD are not decreasing", prefix) 183 } 184 // Check that singular values are non-negative. 185 if minmn > 0 && floats.Min(sAll) < 0 { 186 t.Errorf("Case %v: some singular values from full SVD are negative", prefix) 187 } 188 189 // Do partial SVD and compare the results to sAll, uAll, and vtAll. 190 for _, jobU := range []lapack.SVDJob{lapack.SVDAll, lapack.SVDStore, lapack.SVDOverwrite, lapack.SVDNone} { 191 for _, jobVT := range []lapack.SVDJob{lapack.SVDAll, lapack.SVDStore, lapack.SVDOverwrite, lapack.SVDNone} { 192 if jobU == lapack.SVDOverwrite || jobVT == lapack.SVDOverwrite { 193 // Not implemented. 194 continue 195 } 196 if jobU == lapack.SVDAll && jobVT == lapack.SVDAll { 197 // Already checked above. 198 continue 199 } 200 201 prefix := prefix + ",job=" + svdJobString(jobU) + "U-" + svdJobString(jobVT) + "VT" 202 203 // Restore A to its original values. 204 copy(a, aCopy) 205 206 // Allocate slices for the results of partial SVD and fill them. 207 u := make([]float64, m*ldu) 208 for i := range u { 209 u[i] = rnd.NormFloat64() 210 } 211 vt := make([]float64, n*ldvt) 212 for i := range vt { 213 vt[i] = rnd.NormFloat64() 214 } 215 s := make([]float64, min(m, n)) 216 for i := range s { 217 s[i] = math.NaN() 218 } 219 220 for i := range work { 221 work[i] = math.NaN() 222 } 223 224 ok := impl.Dgesvd(jobU, jobVT, m, n, a, lda, s, u, ldu, vt, ldvt, work, len(work)) 225 if !ok { 226 t.Fatalf("Case %v: unexpected failure in partial Dgesvd", prefix) 227 } 228 229 if minmn == 0 { 230 // No panic and the result is ok, there is 231 // nothing else to check. 232 continue 233 } 234 235 // Check that U has orthogonal columns and that it matches UAll. 236 switch jobU { 237 case lapack.SVDStore: 238 if !hasOrthonormalColumns(blas64.General{Rows: m, Cols: minmn, Data: u, Stride: ldu}) { 239 t.Errorf("Case %v: columns of U are not orthogonal", prefix) 240 } 241 if res := svdPartialUResidual(m, minmn, u, uAll, ldu); res > tol { 242 t.Errorf("Case %v: columns of U do not match UAll", prefix) 243 } 244 case lapack.SVDAll: 245 if !hasOrthonormalColumns(blas64.General{Rows: m, Cols: m, Data: u, Stride: ldu}) { 246 t.Errorf("Case %v: columns of U are not orthogonal", prefix) 247 } 248 if res := svdPartialUResidual(m, m, u, uAll, ldu); res > tol { 249 t.Errorf("Case %v: columns of U do not match UAll", prefix) 250 } 251 } 252 // Check that VT has orthogonal rows and that it matches VTAll. 253 switch jobVT { 254 case lapack.SVDStore: 255 if !hasOrthonormalRows(blas64.General{Rows: minmn, Cols: n, Data: vtAll, Stride: ldvt}) { 256 t.Errorf("Case %v: rows of VT are not orthogonal", prefix) 257 } 258 if res := svdPartialVTResidual(minmn, n, vt, vtAll, ldvt); res > tol { 259 t.Errorf("Case %v: rows of VT do not match VTAll", prefix) 260 } 261 case lapack.SVDAll: 262 if !hasOrthonormalRows(blas64.General{Rows: n, Cols: n, Data: vtAll, Stride: ldvt}) { 263 t.Errorf("Case %v: rows of VT are not orthogonal", prefix) 264 } 265 if res := svdPartialVTResidual(n, n, vt, vtAll, ldvt); res > tol { 266 t.Errorf("Case %v: rows of VT do not match VTAll", prefix) 267 } 268 } 269 // Check that singular values are decreasing. 270 if !sort.IsSorted(sort.Reverse(sort.Float64Slice(s))) { 271 t.Errorf("Case %v: singular values from full SVD are not decreasing", prefix) 272 } 273 // Check that singular values are non-negative. 274 if floats.Min(s) < 0 { 275 t.Errorf("Case %v: some singular values from full SVD are negative", prefix) 276 } 277 if !floats.EqualApprox(s, sAll, tol/10) { 278 t.Errorf("Case %v: singular values differ between full and partial SVD\n%v\n%v", prefix, s, sAll) 279 } 280 } 281 } 282 } 283} 284 285// svdFullResidual returns 286// |A - U*D*VT| / (n * aNorm) 287// where U, D, and VT are as computed by Dgesvd with jobU = jobVT = lapack.SVDAll. 288func svdFullResidual(m, n int, aNorm float64, a []float64, lda int, u []float64, ldu int, d []float64, vt []float64, ldvt int) float64 { 289 // The implementation follows TESTING/dbdt01.f from the reference. 290 291 minmn := min(m, n) 292 if minmn == 0 { 293 return 0 294 } 295 296 // j-th column of A - U*D*VT. 297 aMinusUDVT := make([]float64, m) 298 // D times the j-th column of VT. 299 dvt := make([]float64, minmn) 300 // Compute the residual |A - U*D*VT| one column at a time. 301 var resid float64 302 for j := 0; j < n; j++ { 303 // Copy j-th column of A to aj. 304 blas64.Copy(blas64.Vector{N: m, Data: a[j:], Inc: lda}, blas64.Vector{N: m, Data: aMinusUDVT, Inc: 1}) 305 // Multiply D times j-th column of VT. 306 for i := 0; i < minmn; i++ { 307 dvt[i] = d[i] * vt[i*ldvt+j] 308 } 309 // Compute the j-th column of A - U*D*VT. 310 blas64.Gemv(blas.NoTrans, 311 -1, blas64.General{Rows: m, Cols: minmn, Data: u, Stride: ldu}, blas64.Vector{N: minmn, Data: dvt, Inc: 1}, 312 1, blas64.Vector{N: m, Data: aMinusUDVT, Inc: 1}) 313 resid = math.Max(resid, blas64.Asum(blas64.Vector{N: m, Data: aMinusUDVT, Inc: 1})) 314 } 315 if aNorm == 0 { 316 if resid != 0 { 317 // Original matrix A is zero but the residual is non-zero, 318 // return infinity. 319 return math.Inf(1) 320 } 321 // Original matrix A is zero, residual is zero, return 0. 322 return 0 323 } 324 // Original matrix A is non-zero. 325 if aNorm >= resid { 326 resid = resid / aNorm / float64(n) 327 } else { 328 if aNorm < 1 { 329 resid = math.Min(resid, float64(n)*aNorm) / aNorm / float64(n) 330 } else { 331 resid = math.Min(resid/aNorm, float64(n)) / float64(n) 332 } 333 } 334 return resid 335} 336 337// svdPartialUResidual compares U and URef to see if their columns span the same 338// spaces. It returns the maximum over columns of 339// |URef(i) - S*U(i)| 340// where URef(i) and U(i) are the i-th columns of URef and U, respectively, and 341// S is ±1 chosen to minimize the expression. 342func svdPartialUResidual(m, n int, u, uRef []float64, ldu int) float64 { 343 var res float64 344 for j := 0; j < n; j++ { 345 imax := blas64.Iamax(blas64.Vector{N: m, Data: uRef[j:], Inc: ldu}) 346 s := math.Copysign(1, uRef[imax*ldu+j]) * math.Copysign(1, u[imax*ldu+j]) 347 for i := 0; i < m; i++ { 348 diff := math.Abs(uRef[i*ldu+j] - s*u[i*ldu+j]) 349 res = math.Max(res, diff) 350 } 351 } 352 return res 353} 354 355// svdPartialVTResidual compares VT and VTRef to see if their rows span the same 356// spaces. It returns the maximum over rows of 357// |VTRef(i) - S*VT(i)| 358// where VTRef(i) and VT(i) are the i-th columns of VTRef and VT, respectively, and 359// S is ±1 chosen to minimize the expression. 360func svdPartialVTResidual(m, n int, vt, vtRef []float64, ldvt int) float64 { 361 var res float64 362 for i := 0; i < m; i++ { 363 jmax := blas64.Iamax(blas64.Vector{N: n, Data: vtRef[i*ldvt:], Inc: 1}) 364 s := math.Copysign(1, vtRef[i*ldvt+jmax]) * math.Copysign(1, vt[i*ldvt+jmax]) 365 for j := 0; j < n; j++ { 366 diff := math.Abs(vtRef[i*ldvt+j] - s*vt[i*ldvt+j]) 367 res = math.Max(res, diff) 368 } 369 } 370 return res 371} 372