![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to Dgl.e CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 40BSD / cmd / efl / efltest|
Return to Dgl.e CVS log | Up to [CSRG BSD Unix] / 40BSD / cmd / efl / efltest |
1.1 root 1: Procedure DglsBP(Nu,Order,BC,E)
2:
3: # To determine which ODE should use which boundary condition.
4:
5: # Mnemonic - Double precision Galerkin's method for Linear Systems,
6: # Boundary condition Placement.
7:
8: # Scratch Space Allocated -
9:
10: # S(DglsBP) <= Nu*(4*Nu+15)
11:
12: # Integer words.
13:
14: Integer Nu,Order(Nu,Nu,2),BC(Nu,2,2),E(Nu,2,2)
15:
16: Integer inow,inowold,i,j,l,iMaxord,iCE,Istkgt,iPPS
17: Logical AllZero
18:
19: Struct Nodei { Integer bp,N,j,R(1) }
20:
21: Define Push +1
22: Define Search 0
23: Define Pop -1
24:
25: # Define Node = Is(inow) -> Nodei
26:
27: Include dstack
28:
29: # Check the input for errors.
30:
31: If ( Nu < 1 ) { Seterr("DglsBP - Nu.lt.1",16,1,2) }
32:
33: Do l = 1 , 2
34: {
35:
36: Do i = 1 , Nu
37: {
38: AllZero = .True. # Is Order(i,.,l) = (-1,...,-1)?
39: Do j = 1 , Nu
40: {
41: AllZero &= Order(i,j,l) == (-1)
42: If ( Order(i,j,l) < (-1) | Order(i,j,l) > 2 )
43: { Seterr("DglsBP - Order(i,j,l) not one of -1,0,1,2",41,2,2) }
44: }
45: If ( ( BC(i,1,l) ~= (-2) & BC(i,1,l) ~= 0 ) |
46: ( BC(i,2,l) ~= (-2) & BC(i,2,l) ~= 1 ) )
47: { Seterr("DglsBP - BC(i,.,l) not one of -2,0,1",36,3,2) }
48: If ( AllZero )
49: { Seterr("DglsBP - Order(i,.,l)=(-1,...,-1)",33,4,2) }
50: }
51: }
52:
53: Enter(1)
54:
55: iCE = Istkgt(Nu,2) # Complement of E.
56:
57: # Maxord(i,l) = Max over j=1,...,Nu Order(i,j,l).
58:
59: iMaxord = Istkgt(2*Nu,2)
60: Seti(2*Nu,-1,Is(iMaxord))
61:
62: Do l = 1 , 2
63: {
64: Do i = 1 , Nu
65: {
66: Do j = 1 , Nu
67: {
68: Is(iMaxord+i-1+(l-1)*Nu) = Max0(Is(iMaxord+i-1+(l-1)*Nu),
69: Order(i,j,l))
70: }
71: }
72: }
73:
74: i = 0 ; iPPS = Push
75:
76: While ( i < 4*Nu | iPPS ~= Push )
77: {
78:
79: Switch ( iPPS )
80: {
81:
82: Case Push: # Make a node.
83:
84: inowold = inow
85: i += 1 ; inow = Istkgt(Nu+3,2)
86: Is(inow) -> Nodei.bp = inowold
87:
88: # Get the candidates for E(i).
89:
90: D6lsBP(i,Nu,Order,BC,E,
91: Is(iMaxord),Is(iCE),Is(inow) -> Nodei.R,Is(inow) -> Nodei.N)
92:
93: Is(inow) -> Nodei.j = 0
94: iPPS = Search ; Break
95:
96: Case Search: # Searching a node.
97:
98: Is(inow) -> Nodei.j += 1
99:
100: If ( Is(inow) -> Nodei.j > Is(inow) -> Nodei.N ) # Back-up.
101: { iPPS = Pop ; Next }
102:
103: E(i,1,1) = Is(inow) -> Nodei.R(Is(inow) -> Nodei.j)
104: iPPS = Push ; Break
105:
106: Case Pop: # Backing up a Node.
107:
108: inow = Is(inow) -> Nodei.bp ; Istkrl(1) ; i -= 1
109: iPPS = Search ; Break
110:
111: } # End Switch.
112:
113: If ( i == 0 )
114: { Seterr("DglsBP - Improper Boundary Conditions",37,5,1) ; Break }
115:
116: } # End While.
117:
118: Leave()
119:
120: Return
121:
122: End
123: Procedure D6lsBP(i,Nu,Order,BC,E,
124: Maxord,CE,R,N)
125:
126: Integer i,Nu,Order(Nu,Nu,2),BC(1),E(1), # BC(Nu,2,2),E(Nu,2,2),
127: Maxord(Nu,2),CE(Nu),R(Nu),N # E(i-1),R(N).
128:
129: Integer j,LR,DM,NBCs,l,ii
130:
131: If ( BC(i) < 0 ) { N = 1 ; R(N) = 0 ; Return }
132:
133: # LR = 1 for left, LR = 2 for right.
134:
135: LR = 1+(i-1)/(2*Nu)
136:
137: # DM = 1 for Dirichlet, DM = 2 for Mixed boundary conditions.
138:
139: DM = 1+Mod((i-1)/Nu,2)
140:
141: ii = Mod(i,Nu) ; If ( ii == 0 ) { ii = Nu } # B(i) = B(ii,DM,LR).
142:
143: N = 0
144:
145: Do j = 1 , Nu { CE(j) = j } # CE = Complement of E.
146: If ( i <= 2*Nu )
147: {
148: For ( j = 1 , j < i , j += 1 )
149: {
150: If ( BC(j) >= 0 ) { CE(E(j)) = 0 }
151: }
152: }
153: Else
154: {
155: For ( j = 2*Nu+1 , j < i , j += 1 )
156: {
157: If ( BC(j) >= 0 ) { CE(E(j)) = 0 }
158: }
159: }
160:
161: Do j = 1 , Nu
162: {
163: If ( CE(j) == 0 ) { Next }
164:
165: NBCs = 0
166:
167: For ( l = 1 , l < i , l += 1 )
168: {
169: If ( E(l) == j & BC(l) >= 0 ) { NBCs += 1 }
170: }
171:
172: If ( ( DM == 1 & Maxord(j,LR) > BC(i) ) |
173: ( DM == 2 & Order(j,ii,LR) > BC(i) ) )
174: {
175: If ( NBCs < Max0(Maxord(j,1),Maxord(j,2)) )
176: { N += 1 ; R(N) = j }
177: }
178: }
179:
180: Return
181:
182: End
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.