forked from cpmech/gosl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
densesol.go
116 lines (103 loc) · 2.78 KB
/
densesol.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
106
107
108
109
110
111
112
113
114
115
116
// Copyright 2016 The Gosl Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package la
import (
"math"
"github.com/cpmech/gosl/chk"
"github.com/cpmech/gosl/la/oblas"
)
// DenSolve solves dense linear system using LAPACK (OpenBLAS)
//
// Given: A ⋅ x = b find x such that x = A⁻¹ ⋅ b
//
func DenSolve(x Vector, A *Matrix, b Vector, preserveA bool) {
a := A
if preserveA {
a = NewMatrix(A.M, A.N)
copy(a.Data, A.Data)
}
copy(x, b)
ipiv := make([]int32, A.M)
oblas.Dgesv(A.M, 1, a.Data, A.M, ipiv, x, A.M)
}
// Cholesky returns the Cholesky decomposition of a symmetric positive-definite matrix
//
// a = L * trans(L)
//
func Cholesky(L, a *Matrix) {
for j := 0; j < a.M; j++ { // loop over columns
for i := j; i < a.M; i++ { // loop over lower diagonal rows (including diagonal)
amsum := a.Get(i, j)
for k := 0; k < j; k++ {
amsum -= L.Get(i, k) * L.Get(j, k)
}
if i == j {
if amsum <= 0.0 {
chk.Panic("Cholesky factorization failed due to non positive-definite matrix")
}
L.Set(i, j, math.Sqrt(amsum))
} else {
L.Set(i, j, amsum/L.Get(j, j))
}
}
}
}
// SolveRealLinSysSPD solves a linear system with real numbres and a Symmetric-Positive-Definite (SPD) matrix
//
// x := inv(a) * b
//
// NOTE: this function uses Cholesky decomposition and should be used for small systems
func SolveRealLinSysSPD(x Vector, a *Matrix, b Vector) {
// Cholesky factorisation
L := NewMatrix(a.M, a.M)
Cholesky(L, a)
// solve L*y = b storing y in x
for i := 0; i < a.M; i++ {
bmsum := b[i]
for k := 0; k < i; k++ {
bmsum -= L.Get(i, k) * x[k]
}
x[i] = bmsum / L.Get(i, i)
}
// solve trans(L)*x = y with y==x
for i := a.M - 1; i >= 0; i-- {
bmsum := x[i]
for k := i + 1; k < a.M; k++ {
bmsum -= L.Get(k, i) * x[k]
}
x[i] = bmsum / L.Get(i, i)
}
}
// SolveTwoRealLinSysSPD solves two linear systems with real numbres and Symmetric-Positive-Definite (SPD) matrices
//
// x := inv(a) * b and X := inv(a) * B
//
// NOTE: this function uses Cholesky decomposition and should be used for small systems
func SolveTwoRealLinSysSPD(x, X Vector, a *Matrix, b, B Vector) {
// Cholesky factorisation
L := NewMatrix(a.M, a.M)
Cholesky(L, a)
// solve L*y = b storing y in x
for i := 0; i < a.M; i++ {
bmsum := b[i]
Bmsum := B[i]
for k := 0; k < i; k++ {
bmsum -= L.Get(i, k) * x[k]
Bmsum -= L.Get(i, k) * X[k]
}
x[i] = bmsum / L.Get(i, i)
X[i] = Bmsum / L.Get(i, i)
}
// solve trans(L)*x = y with y==x
for i := a.M - 1; i >= 0; i-- {
bmsum := x[i]
Bmsum := X[i]
for k := i + 1; k < a.M; k++ {
bmsum -= L.Get(k, i) * x[k]
Bmsum -= L.Get(k, i) * X[k]
}
x[i] = bmsum / L.Get(i, i)
X[i] = Bmsum / L.Get(i, i)
}
}