/
operations.go
125 lines (112 loc) · 2.27 KB
/
operations.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
117
118
119
120
121
122
123
124
125
package sparse
import (
"fmt"
"github.com/Konstantin8105/errors"
)
// IsSym return true if matrix A is symmetrical
func IsSym(A *Matrix) (ok bool, err error) {
// check input data
et := errors.New("Function Limits: check input data")
if A == nil {
_ = et.Add(fmt.Errorf("matrix A is nil"))
}
if A != nil && A.nz != -1 {
_ = et.Add(fmt.Errorf("matrix A is not CSC(Compressed Sparse Column) format"))
}
if A != nil {
if A.n == 0 && A.m == 0 {
_ = et.Add(fmt.Errorf("matrix A is empty"))
}
if A.n != A.m {
_ = et.Add(fmt.Errorf("matrix A is not square"))
}
}
if et.IsError() {
return false, et
}
// initialization
n, Ap, Ai, Ax := A.n, A.p, A.i, A.x
// check amount upper and lower entries
{
nzUpper := 0
nzLower := 0
for j := 0; j < n; j++ {
for p := Ap[j]; p < Ap[j+1]; p++ {
i := Ai[p] // row
j := j // column
if i == j {
// ignore diagonal
continue
}
if i > j { // lower
nzLower++
continue
}
// upper
nzUpper++
}
}
if nzUpper != nzLower {
return false, fmt.Errorf(
"amount entries on up and low diagonal is not same: %d != %d",
nzUpper, nzLower)
}
}
defer func() {
// unmark all
for j := 0; j < n; j++ {
for p := Ap[j]; p < Ap[j+1]; p++ {
Ai[p] = unflip(Ai[p])
}
}
}()
// check values
for j := 0; j < n; j++ {
for p := Ap[j]; p < Ap[j+1]; p++ {
if marked(Ai, p) {
// ignore marked
continue
}
i := Ai[p] // row
j := j // column
x := Ax[p] // value
if i == j {
continue
}
// finding
var found bool
for p2 := Ap[i]; p2 < Ap[i+1]; p2++ {
i2 := Ai[p2] // row
j2 := i // column
x2 := Ax[p2] // value
if marked(Ai, p2) || i != j2 || j != i2 || i2 == j2 {
// ignore marked
// or
// coordinates of entries is not same
// or
// ignore diagonal entries
continue
}
found = true
if x != x2 {
return false, fmt.Errorf(
"matrix is not symmetric. Value[%d,%d] = %.5e. Value[%d,%d] = %5e.",
i, j, x,
i2, j2, x2)
}
if found {
// mark
mark(Ai, p2)
mark(Ai, p)
break
}
}
if !found {
return false, fmt.Errorf(
"matrix is not symmetric. Cannot find entry [%d,%d]",
i, j)
}
}
}
return true, nil
}