forked from hrautila/linalg
/
sytrf.go
105 lines (93 loc) · 2.79 KB
/
sytrf.go
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
// Copyright (c) Harri Rautila, 2012,2013
// This file is part of github.com/hrautila/linalg/lapack package.
// It is free software, distributed under the terms of GNU Lesser General Public
// License Version 3, or any later version. See the COPYING tile included in this archive.
package lapack
import (
//"errors"
"fmt"
"github.com/hrautila/linalg"
"github.com/hrautila/matrix"
)
/*
LDL^T factorization of a real or complex symmetric matrix.
PURPOSE
Computes the LDL^T factorization of a real or complex symmetric
n by n matrix A. On exit, A and ipiv contain the details of the
factorization.
ARGUMENTS
A float or complex matrix
ipiv int vector of length at least n
OPTIONS
uplo PLower or PUpper
n nonnegative integer. If negative, the default value is used.
ldA positive integer. ldA >= max(1,n). If zero, the default
value is used.
offsetA nonnegative integer;
*/
func Sytrf(A matrix.Matrix, ipiv []int32, opts ...linalg.Option) error {
switch A.(type) {
case *matrix.FloatMatrix:
return SytrfFloat(A.(*matrix.FloatMatrix), ipiv, opts...)
case *matrix.ComplexMatrix:
return SytrfComplex(A.(*matrix.ComplexMatrix), ipiv, opts...)
}
return onError("Sytrf: unknown types")
}
func SytrfFloat(A *matrix.FloatMatrix, ipiv []int32, opts ...linalg.Option) error {
pars, err := linalg.GetParameters(opts...)
if err != nil {
return err
}
ind := linalg.GetIndexOpts(opts...)
err = checkSytrf(ind, A, ipiv)
if err != nil {
return err
}
if ind.N == 0 {
return nil
}
Aa := A.FloatArray()
uplo := linalg.ParamString(pars.Uplo)
info := dsytrf(uplo, ind.N, Aa[ind.OffsetA:], ind.LDa, ipiv)
if info != 0 {
return onError(fmt.Sprintf("Sytrf: lapack error %d", info))
}
return nil
}
func SytrfComplex(A *matrix.ComplexMatrix, ipiv []int32, opts ...linalg.Option) error {
return onError(fmt.Sprintf("Sytrf: complex not yet implemented"))
}
func checkSytrf(ind *linalg.IndexOpts, A matrix.Matrix, ipiv []int32) error {
arows := ind.LDa
if ind.N < 0 {
ind.N = A.Rows()
if ind.N != A.Cols() {
return onError("A not square")
}
}
if ind.N == 0 {
return nil
}
if ind.LDa == 0 {
ind.LDa = max(1, A.LeadingIndex())
arows = max(1, A.Rows())
}
if ind.LDa < max(1, ind.N) {
return onError("Sytrf: lda")
}
if ind.OffsetA < 0 {
return onError("Sytrf: offsetA")
}
sizeA := A.NumElements()
if sizeA < ind.OffsetA+(ind.N-1)*arows+ind.N {
return onError("Sytrf: sizeA")
}
if ipiv != nil && len(ipiv) < ind.N {
return onError("Sytrf: size ipiv")
}
return nil
}
// Local Variables:
// tab-width: 4
// End: