-
Notifications
You must be signed in to change notification settings - Fork 1
/
band.m
95 lines (87 loc) · 1.76 KB
/
band.m
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
function [C,E] = band(J,A,B,C,D,E,G,X,Y,N,NJ)
NP1 = N+1;
if J == 1
for I = 1:N
D(I,2*N+1) = G(I);
for L = 1:N
LPN = L+N;
D(I,LPN) = X(I,L);
end
end
[det,D] = matinv(N,2*N+1,B,D);
if det == 0
fprintf('Determinant = 0 at j = %d\n',J)
return
else
for K = 1:N
E(K,NP1,1) = D(K,2*N+1);
for L = 1:N
E(K,L,1) = -D(K,L);
LPN = L+N;
X(K,L) = -D(K,LPN);
end
end
return
end
end
if J == 2
for I = 1:N
for K = 1:N
for L = 1:N
D(I,K) = D(I,K) + A(I,L)*X(L,K);
end
end
end
end
if J == NJ
for I = 1:N
for L = 1:N
G(I) = G(I)-Y(I,L)*E(L,NP1,J-2);
for M = 1:N
A(I,L) = A(I,L)+Y(I,M)*E(M,L,J-2);
end
end
end
end
for I = 1:N
D(I,NP1) = -G(I);
for L = 1:N
D(I,NP1) = D(I,NP1)+A(I,L)*E(L,NP1,J-1);
for K = 1:N
B(I,K) = B(I,K)+A(I,L)*E(L,K,J-1);
end
end
end
[det,D] = matinv(N,NP1,B,D);
if det == 0
fprintf('Determinant = 0 at j = %d\n',J)
return
else
for K = 1:N
for M = 1:NP1
E(K,M,J) = -D(K,M);
end
end
end
if J < NJ
return
else
for K = 1:N
C(K,J) = E(K,NP1,J);
end
for JJ = 2:NJ
M = NJ-JJ+1;
for K = 1:N
C(K,M) = E(K,NP1,M);
for L = 1:N
C(K,M) = C(K,M)+E(K,L,M)*C(L,M+1);
end
end
end
for L = 1:N
for K = 1:N
C(K,1) = C(K,1)+X(K,L)*C(L,3);
end
end
end
end