More matrix fun

[david_uwi]

As we seem to be looking at matrix manipulation...
Here's a program I wrote nearly 30 years ago. It diagonalizes a real symmetric matrix. This gives (in maths-speak) the eigenvalues and eigenvectors. It''s only set up for a 3x3 (which is read from a data file in the order 1,1 2,2 3,3 1,2 1,3 2,3),
but it will do any size (but note the time taken is approx N cubed where N is the size of the matrix). Enjoy...

Code: QB64: [Select]

1. 5 DEFINT I-N
2. n2 = 3
3. DIM e(n2, n2), X(n2, n2)
4. OPEN "c:\cw1\mymatrix.txt" FOR INPUT AS #1
5. INPUT #1, e(1, 1), e(2, 2), e(3, 3), e(1, 2), e(1, 3), e(2, 3)
6. CLOSE #1
7. e(2, 1) = e(1, 2): e(3, 2) = e(2, 3): e(3, 1) = e(1, 3)
8. PRINT "input matrix"
9. FOR i = 1 TO n2
10.     FOR J = 1 TO n2
11.         PRINT USING "####.#####"; e(i, J);
12.     NEXT J
13.     PRINT
14. NEXT i
15. PRINT
16.  
17. GOSUB 4000
18. PRINT "EIGENVALUES"
19. FOR i = 1 TO n2
20.     FOR J = 1 TO n2
21.         PRINT USING "####.#####"; e(i, J);
22.     NEXT J
23.     PRINT
24. NEXT i
25. PRINT
26. PRINT "EIGENVECTORS"
27. FOR i = 1 TO n2
28.     FOR J = 1 TO n2
29.         PRINT USING "####.#####"; X(i, J);
30.     NEXT J
31.     PRINT
32. NEXT i
33.  
34. CLOSE #1
35. GOTO 130
36. 100 PRINT "*********************"
37. 110 PRINT "INCORRECT INPUT!!!!!"
38. 120 PRINT "*********************"
39. 130 END
40. 4000 REM DIAGONALISATION USING JACOBI'S METHOD (D.S. 1993)
41. 4005 M = (n2 * n2 - n2) / 2
42. 4010 TRA = M * .00001
43. FOR i = 1 TO n2
44.     FOR J = 1 TO n2
45.         X(i, J) = 0!
46.     NEXT J
47.     X(i, i) = 1!
48. NEXT i
49. 4090 FOR L = 1 TO M * 10
50.    4100 Q = 0: AM = -1
51.    4110 FOR J = 1 TO n2 - 1
52.        4120 FOR i = J + 1 TO n2
53.            4130 Z = ABS(e(i, J))
54.            4140 Q = Q + Z
55.            4150 IF Z < AM GOTO 4170
56.            4160 AM = Z: IX = i: JX = J
57.        4170 NEXT i
58.    4180 NEXT J
59.    4190 IF Q < TRA THEN RETURN
60.    4200 A1 = e(IX, IX) - e(JX, JX)
61.    4210 A5 = e(IX, JX)
62.    4220 IF ABS(A5) * 50 > ABS(A1) GOTO 4270
63.    4230 A6 = 2 * A1 * A1 + 3 * A5 * A5
64.    4240 C = -2 * A1 * A5 / A6
65.    4250 s = 1 - A5 * A5 / A6
66.    4260 GOTO 4320
67.    4270 A2 = SQR(A1 * A1 + 4 * A5 * A5)
68.    4280 IF A1 > 0 THEN A2 = -A2
69.    4290 A3 = (A1 + A2) / (2 * A5)
70.    4300 s = 1 / SQR(A3 * A3 + 1)
71.    4310 C = A3 * s
72.    4320 FOR J = 1 TO n2
73.        4330 TM = e(JX, J) * s + e(IX, J) * C
74.        4340 e(IX, J) = e(IX, J) * s - e(JX, J) * C
75.        4350 e(JX, J) = TM
76.        4360 TM = X(JX, J) * s + X(IX, J) * C
77.        4370 X(IX, J) = X(IX, J) * s - X(JX, J) * C
78.        4380 X(JX, J) = TM
79.    4390 NEXT J
80.    4400 FOR J = 1 TO n2
81.        4410 TM = e(J, JX) * s + e(J, IX) * C
82.        4420 e(J, IX) = e(J, IX) * s - e(J, JX) * C
83.        4430 e(J, JX) = TM
84.    4440 NEXT J
85. 4450 NEXT L
86. 4460 PRINT "CONVERGENCE FAILURE": CLOSE #1
87. 4490 END

[david_uwi]

That's amazing the Jacobi diagonalization is very similar to mine. Maybe I was influenced by "Numerical Methods" though some of the programs in the book are very poor efforts.
I would claim that mine is better as it avoids rounding errors by moving from using a square root solution as the off-diagonals get small (notice I only use single precision). There is also a problem with the solution "stagnating" and oscillating around a solution which mine overcomes. I often diagonalize several thousand matrices and I can't afford for some to not converge or converge very slowly. Oh and mine has the bonus of line numbers and GOTOs.
I'm looking at the "Hermitian routine" (that does matrices with a+ib on one side and a-ib on the other side of the diagonal).
The trick is that the result (eigenvalues) has to be real so you can initially reduce the matrix to a real symmetric and go from there.

[david_uwi]

If you want a real world example of this routine in action. Here is my NMR (nuclear magnetic resonance) spectrum calculator. In NMR you are primarily looking at hydrogen and their resonance frequency depends on the surrounding atoms (this is called chemical shift) there is also interaction between the hydrogens (so called spin-spin coupling). It is this spin-spin coupling that produces off-diagonal elements in the matrix -that need to be removed. I have given a small demo data file for 1,2-dichlorobenzene a very simple molecule that surprisingly produces a rather complex spectrum.
This is the data file (if you look at the program this is expected in c:\cw1\o-dichlorobenzene.txt). Almost forgot the program is also expecting to be able to load fonts from c:\windows\fonts.

dichlorobenzene
4,100
6.0,6.5,6.5,6.0
5,1,0,5,1,5
5.5,7.0,.5,0

line one is just the title
line 2 is the number of hydrogens(4), and the spectrometer operating frequency (100 MHz)
line 3 is the chemical shifts (4 hydrogens 4 shifts, but 1 and 4 are the same as are 2 and 3)
then come the couplings in the order 1,2 1,3 1,4 2,3 2,4 3,4
line 4 is for the output low chemical shift, high, line broadening factor, integration (which is set off)

Code: QB64: [Select]

1. DEFINT I-N
2. DEFDBL A-H, O-Z
3. OPEN "I", #1, "c:\cw1\o-dichlorobenzene.txt"
4. INPUT #1, b$
5. INPUT #1, n, spop
6. IM = 2 ^ n
7. DIM AH(IM, n) AS INTEGER, AX(IM, IM) AS INTEGER
8. DIM TP(4096), F(4096), X(IM, IM), E(IM, IM)
9. DIM aj(28, 28), VH(n), XX(4000), Y(4000), YD(30), XD(40)
10. LINE INPUT #1, Z$ 'INPUT OF CHEMICAL SHIFTS
11. GOSUB charsep1
12. REM INPUT OF SPIN-SPIN COUPLING J12,J13...J23,J24..etc
13. LINE INPUT #1, Z$
14. GOSUB charsep2
15. REM freq min,freq max, line width,integration (0/1)
16. INPUT #1, fmin, fmax, t2, itg
17. fmin = fmin * spop: fmax = fmax * spop
18. TK = TIMER
19. REM CALCULATING ELEMENTS IN HAMILTONIAN MATRIX
20. xsca = 960: ysca = 480 'spectrum size
21. N2 = IM
22. FOR i = 0 TO N2 - 1
23.     II = i
24.     FOR IN = n - 1 TO 0 STEP -1
25.         JJ = 2 ^ IN
26.         IF (II - JJ) < 0 THEN TJ = -1 ELSE TJ = 1
27.         IF TJ = 1 THEN II = II - JJ
28.         AH(i + 1, IN + 1) = TJ
29.     NEXT IN
30. NEXT i
31. FOR i = 1 TO N2
32.     E(i, i) = 0: AX(i, i) = 0
33.     FOR L = 1 TO n
34.         E(i, i) = E(i, i) + VH(L) * .5 * AH(i, L)
35.         TD = AH(i, L) * .25
36.         FOR k = L + 1 TO n
37.             E(i, i) = E(i, i) + aj(L, k) * TD * AH(i, k)
38.         NEXT k
39.     NEXT L
40. NEXT i
41. FOR i = 1 TO N2
42.     X(i, i) = 1
43.     FOR L = i + 1 TO N2
44.         MA = 0: MB = 0: K1 = 0: K2 = 0
45.         X(i, L) = 0: X(L, i) = 0
46.         FOR k = 1 TO n
47.             TG = AH(i, k) * AH(L, k) - 1
48.             MA = MA + TG
49.             MB = MB + TG * AH(L, k)
50.             IF TG <> 0 THEN
51.                 IF K1 > 0 THEN K2 = k ELSE K1 = k
52.             END IF
53.         NEXT k
54.         IF MA = -4 AND MB = 0 THEN
55.             E(i, L) = .5 * aj(K1, K2)
56.             E(L, i) = E(i, L)
57.         END IF
58.     NEXT L
59. NEXT i
60. GOSUB jacobi
61. 'E matrix has the eigenvalues (energies)
62. 'X matrix has the eigenvectors (wavefunction mixtures)
63. IAY = 0
64. FOR i = 1 TO N2
65.     FOR j = i + 1 TO N2
66.         FOR k = 1 TO n
67.             IAY = IAY + AH(i, k) * AH(j, k)
68.         NEXT k
69.         IF IAY = n - 2 THEN AX(i, j) = 1 ELSE AX(i, j) = 0
70.         AX(j, i) = AX(i, j)
71.         IAY = 0
72.     NEXT j
73. NEXT i
74. m = 0
75. FOR L = 1 TO N2 - 1 'we have 4 nested loops - this is slow!
76.     FOR k = L + 1 TO N2
77.         ff = ABS(E(k, k) - E(L, L)) 'frequencies
78.         IF ABS(ff) > 50000 GOTO 10
79.         tpp = 0
80.         FOR j = 1 TO N2
81.             IF ABS(X(L, j)) > .01 THEN
82.                 FOR i = 1 TO N2
83.                     IF AX(j, i) <> 0 THEN tpp = tpp + X(k, i) * X(L, j)
84.                 NEXT i
85.             END IF
86.         NEXT j
87.         IF ABS(tpp) > .02 THEN
88.             m = m + 1
89.             F(m) = ff
90.             TP(m) = tpp * tpp 'transition probabilities
91.         END IF
92.    10 NEXT k
93. NEXT L
94. CLOSE #1
95. GOSUB sort
96. GOSUB lorentz
97. GOSUB outscr
98. CLOSE #1
99. END 'end of main
100. sort:
101. FOR i = 1 TO m - 1
102.     FM = F(i): JM = i
103.     FOR j = i + 1 TO m
104.         IF F(j) > FM THEN FM = F(j): JM = j
105.     NEXT j
106.     F(JM) = F(i): F(i) = FM
107.     TPT = TP(JM)
108.     TP(JM) = TP(i): TP(i) = TPT
109. NEXT i
110. RETURN
111. lorentz:
112. REM CONSTRUCTING LORENTZIANS
113. 'lineshapes are lorentzian
114. SW = fmax - fmin
115. MM = 3600
116. AI = SW / MM
117. FOR N5 = 1 TO MM + 1
118.     XX(N5) = fmin + AI * (N5 - 1)
119.     Y(N5) = 0
120. NEXT N5
121. WID = 4 / t2
122. TA = 4 * 3.14159 * 3.14159 * t2 * t2
123. TB = WID / AI
124. FOR N5 = 1 TO m
125.     C2 = (F(N5) - fmin) / AI
126.     IF ABS(F(N5)) > 10000 GOTO 100 'this really needs to include spop
127.     IF ABS(C2 - TB) > 30000 THEN IJC = 30000
128.     IF IJC <> 30000 THEN IJC = INT(C2 - TB)
129.     IF IJC < 1 THEN IJC = 1
130.     IF IJC < MM THEN
131.         IJD = INT(C2 + TB)
132.         IF IJD > MM THEN IJD = MM
133.         C1 = 2 * t2 * TP(N5)
134.         FOR MI = IJC TO IJD
135.             C3 = F(N5) - XX(MI)
136.             Y(MI) = Y(MI) + C1 / (1 + TA * C3 * C3)
137.         NEXT MI
138.     END IF
139. 100 NEXT N5
140. RETURN
141. jacobi:
142. REM DIAGONALISATION USING JACOBI'S METHOD (D.S.)
143. 'removes biggest off-diagonal elements first
144. m = (N2 * N2 - N2) / 2
145. TRA = m * .001
146. FOR L = 1 TO m * 2
147.     Q = 0: AM = -1
148.     FOR j = 1 TO N2 - 1
149.         FOR i = j + 1 TO N2
150.             Z = ABS(E(i, j))
151.             Q = Q + Z
152.             IF Z > AM THEN AM = Z: IX = i: JX = j
153.         NEXT i
154.     NEXT j
155.     IF Q < TRA THEN RETURN
156.     A1 = E(IX, IX) - E(JX, JX)
157.     A5 = E(IX, JX)
158.     IF ABS(A5) * 50 < ABS(A1) THEN
159.         A6 = 2 * A1 * A1 + 3 * A5 * A5
160.         C = -2 * A1 * A5 / A6 'this is not the best approx.
161.         s = 1 - A5 * A5 / A6 'but it stops the routine stagnating
162.     ELSE
163.         A2 = SQR(A1 * A1 + 4 * A5 * A5)
164.         IF A1 > 0 THEN A2 = -A2
165.         A3 = (A1 + A2) / (2 * A5)
166.         s = 1 / SQR(A3 * A3 + 1)
167.         C = A3 * s
168.     END IF
169.     FOR j = 1 TO N2
170.         TM = E(JX, j) * s + E(IX, j) * C
171.         E(IX, j) = E(IX, j) * s - E(JX, j) * C
172.         E(JX, j) = TM
173.         TM = X(JX, j) * s + X(IX, j) * C
174.         X(IX, j) = X(IX, j) * s - X(JX, j) * C
175.         X(JX, j) = TM
176.     NEXT j
177.     FOR j = 1 TO N2
178.         TM = E(j, JX) * s + E(j, IX) * C
179.         E(j, IX) = E(j, IX) * s - E(j, JX) * C
180.         E(j, JX) = TM
181.     NEXT j
182. NEXT L
183. 'program ends up here if the matrix does not
184. 'diagonalize in 2M iterations
185. PRINT "CONVERGENCE FAILURE": CLOSE #1
186. END
187. charsep1:
188. REM SEPARATING INPUT CHARACTER INTO NUMBERS
189. FOR L = 1 TO n - 1
190.     II = LEN(Z$)
191.     JJ = INSTR(Z$, ",")
192.     VH(L) = VAL(LEFT$(Z$, JJ - 1)) * spop
193.     Z$ = RIGHT$(Z$, II - JJ)
194. NEXT L
195. VH(n) = VAL(Z$) * spop
196. RETURN
197. END
198. charsep2:
199. REM SEPARATING INPUT CHARACTER INTO NUMBERS
200. FOR i = 1 TO n - 2
201.     FOR j = i + 1 TO n
202.         II = LEN(Z$)
203.         JJ = INSTR(Z$, ",")
204.         aj(i, j) = VAL(LEFT$(Z$, JJ - 1))
205.         Z$ = RIGHT$(Z$, II - JJ)
206.     NEXT j
207. NEXT i
208. aj(n - 1, n) = VAL(Z$)
209. RETURN
210. END
211. outscr:
212. SCREEN _NEWIMAGE(xsca + 100, ysca + 120, 2)
213. font1& = _LOADFONT("c:\windows\fonts\times.ttf", 32, "bold")
214. _FONT font1&
215.  
216. YMIN = 5000!: YMAX = -1!: xmin = 5000!: xmax = 0!
217. C1 = (fmax - fmin) / 3600: i = 3600
218. YQM = 0
219. FOR k = 1 TO i
220.     IF YMIN > Y(k) THEN YMIN = Y(k)
221.     IF YMAX < Y(k) THEN KMAX = k
222.     IF YMAX < Y(k) THEN YMAX = Y(k)
223.     XX(k) = (C1 * k + fmin) / spop
224.     IF xmin > XX(k) THEN xmin = XX(k)
225.     IF xmax < XX(k) THEN xmax = XX(k)
226.     YQM = YQM + Y(k)
227. NEXT k
228. XO = XX(i): YO = Y(i)
229. XO = ((xmax - XX(i)) / (xmax - xmin) * xsca) + 40
230. YO = ((Y(i) - YMIN) / (YMAX - YMIN) * ysca) - 20
231. REM XO = INT(XO + .5): YO = INT(YO + .5)
232. YYO = 100: YYY = 100
233. FOR k = i - 1 TO 1 STEP -1
234.     XX = ((xmax - XX(k)) / (xmax - xmin) * xsca) + 40
235.     YY = ((Y(k) - YMIN) / (YMAX - YMIN) * ysca) - 20
236.     YN = (YY + YO * TC) / (1 + TC)
237.     REM XX = INT(XX + .5): YY = INT(YN + .5)
238.     LINE (XO, ysca - YO)-(XX, ysca - YY)
239.     YYY = YYY + Y(k) * ysca * .8 / YQM
240.     'integration line
241.     IF itg = 1 THEN LINE (XO, ysca - YYO)-(XX, ysca - YYY)
242.     YYYO = YYY
243.     XO = XX: YO = YN: YYO = YYY
244. NEXT k
245. LINE (40, 10)-(40, ysca + 40)
246. LINE (40, ysca + 40)-(xsca + 40, ysca + 40)
247. LINE (40, 10)-(xsca + 40, 10)
248. LINE (xsca + 40, ysca + 40)-(xsca + 40, 10)
249. kg = -1
250. 20 k = 0
251. kg = kg + 1
252. g(0) = .05: g(1) = .1: g(2) = .2: g(3) = .5: g(4) = 1: g(6) = 20: g(7) = 50: g(8) = 100: g(9) = 200
253. IF xmin = 0 THEN xs = -g(kg)
254. IF xmin <> 0 THEN xs = xmin - g(kg)
255. xi = xs
256. 30 xi = xi + g(kg)
257. 'PRINT xi; g(kg); k; xmin; xmax; xs; kg: INPUT sa$
258. IF xi <= xmax + .001 THEN
259.     k = k + 1
260.     IF k > 10 GOTO 20
261.     XD(k) = xi
262.     GOTO 30
263. END IF
264. XD(k + 1) = xi
265. FOR i = 1 TO k + 1
266.     XXD = (((xmax - XD(i)) / (xmax - xmin)) * xsca) + 40
267.     XC = INT(((XXD - 10) / xsca) * xsca)
268.     XXD = INT(XXD + .5)
269.     cc$ = LTRIM$(STR$((INT(XD(i) * 100 + .5) / 100)))
270.     IF LEN(cc$) = 1 THEN cc$ = cc$ + ".0"
271.     _PRINTSTRING (XC, ysca + 50), cc$
272.     LINE (XXD, ysca + 40)-(XXD, ysca + 50)
273. 40 NEXT i
274. A$ = "Chemical Shift ë (ppm)"
275. _PRINTSTRING (xsca / 2 - 100, ysca + 85), A$
276. ypos = (ytop - ysca) / 8
277. xpos = xsca / 12
278. LOCATE 1, xpos + 10
279. DO
280.     h$ = INKEY$
281. LOOP WHILE h$ = ""
282. IF UCASE$(h$) = "P" THEN GOSUB onebit
283. 'input p to store as a bitmap
284. 7510 RETURN
285.  
286. onebit:
287. '1-bit bitmap loader
288. OPEN b$ + ".bmp" FOR BINARY AS #1
289. pow2%(0) = 1
290. FOR j% = 1 TO 7
291.     pow2%(j%) = pow2%(j% - 1) + pow2%(j% - 1)
292. NEXT j%
293. filetype$ = "BM"
294. offset& = 62&
295. pixwidth& = 1600&
296. pixht& = 820&
297. ppbyte& = 8&
298. header& = 40&
299. datsize& = (pixwidth& * pixht& / ppbyte&)
300. filesize& = datsize& + offset&
301. reserved% = 0
302. planes% = 1
303. bpp% = 1
304. compress& = 0&
305. ppm& = 3880&
306. numcol& = 2&
307. impcol& = 0&
308. 'storing file header
309. PUT #1, , filetype$
310. PUT #1, , filesize&
311. PUT #1, , reserved%
312. PUT #1, , reserved%
313. PUT #1, , offset&
314. 'storing DIB header
315. PUT #1, , header&
316. PUT #1, , pixwidth&
317. PUT #1, , pixht&
318. PUT #1, , planes%
319. PUT #1, , bpp%
320. PUT #1, , compress&
321. PUT #1, , datsize&
322. PUT #1, , ppm&
323. PUT #1, , ppm&
324. PUT #1, , numcol&
325. PUT #1, , impcol&
326. 'storing black and white
327. black$ = CHR$(0) + CHR$(0) + CHR$(0) + CHR$(0)
328. white$ = CHR$(255) + CHR$(255) + CHR$(255) + CHR$(0)
329. PUT #1, , black$
330. PUT #1, , white$
331. 'storing bitmap
332. FOR yvalue% = CINT(pixht&) - 1 TO 0 STEP -1
333.     FOR xvalue% = 0 TO CINT(pixwidth&) - 1 STEP CINT(ppbyte&)
334.         FOR bits% = 0 TO 7
335.             pix% = POINT(xvalue% + bits%, yvalue%)
336.             'change to pix% > 0 to get white on black
337.             IF pix% = 0 THEN pixel% = (pixel% OR pow2%(7 - bits%))
338.         NEXT bits%
339.         qpix$ = CHR$(pixel%)
340.         pixel% = 0
341.         PUT #1, , qpix$
342.     NEXT xvalue%
343. NEXT yvalue%
344. CLOSE #1
345. RETURN
346.  
347.  
348.  

[bplus]

Oh and mine has the bonus of line numbers and GOTOs.

I find it amusing to consider it a bonus. I have to ask, sorry inquiring minds have to know, why is that better?

To me GOTO is completely neutral thing (because I am an enlightened being LOL) and line numbers a nuisance for coding, but maybe I will get insight or find out you are joking around a little :-)

[david_uwi]

I am just amazed at the backward compatibility. All other programming languages there is a good chance you code will stop working following every upgrade.
I used to do a bit of FORTRAN (when I was but a young lad) , but modern FORTRAN is totally different (and to me indecipherable). You would even struggle to get an old quickbasic program to run on freebasic (but not on QB64).
Notice in the second program the GOTOs have gone (and most of the line numbers) -Mac (rip) showed me the error of my ways.

[bplus]

Ah for backwards compatibility, OK yes, saves allot of work updating old working programs.

Can't see anyone wanting to make new programs requiring line numbers.

[NOVARSEG]

Was it GWBASIC that had GOTO line number?

Did GWBASIC have labels?   I can't remember

[david_uwi]

Yes GWBasic had it all, but it was surprisingly slow. I believe the WHILE..WEND is also for back compatibility (when you could easily use a DO..LOOP).

[bplus]

Only line numbers, no line labels (strings)
http://www.antonis.de/qbebooks/gwbasman/index.html