Author |
Topic: Cholesky solution of matrices (Read 542 times) |
|
Himself
New Member
member is offline
Posts: 3
|
 |
Re: Cholesky solution of matrices
« Reply #2 on: Jun 11th, 2017, 1:56pm » |
|
:
REM Cholesky App
:
REM Andre-Louis Cholesky (1885 - 1918) : French military officer and mathematician.
REM Test data extracted from J. C. Nash 1979 publication: who used a HP 9830 workstation.
REM
REM The Cholesky algorithm is here tittivated for precise Engineering Surveying applications,
REM and others busying themselves with multi-parameter least squares calculations. It is
REM space-saving and speedy.
REM
REM The advent of the *FLOAT 80 function enhances the validity of the "most likely" solution
REM for the value of the parameters ( to 19 significant figures no less ).
REM
REM Modern surveying instruments called Total Stations can provide sub-second and sub-millimetre measurement values.
REM This requires computational methods to match ongoing technological developments and engineering challenges,
REM without the arithmetical burden, and limitations, of the matrix inversion approach to the sought-after solution.
REM
REM Run under Windows 10 Pro (1607): BB4W v 6.10a : Computer HP Z400, Quad core 2.40 Xeon.
REM
:
*FLOAT 80
*FONT Courier,16,B
:
INSTALL @lib$+"FNUSING"
:
(Start)
INPUT " For this demonstrationn enter 5 as the number of unknown parameters for this test";n%
IF n% <> 5 THEN
SOUND 1,-15,52,-1
SOUND 1,-15,56,50
WAIT 50
SOUND &11,-15,100,10
CLS
GOTO (Start)
ENDIF
:
: REM Normal equations n% x n% symmetric, positive definite.
:
DIM a(n%*(n%+1)/2) : REM Vector storage (n*(n+1))/2 where n equals the
DIM x(n%) : REM order of the matrix.
DIM q%(n%)
: REM Vector storage of the lower triangle of the normal equations.
a(1)=18926823
a(2)=6359705 : a(3)=2164379
a(4)=10985647 : a(5)=3734131 : a(6)=6445437
a(7)=3344971 : a(8)=1166559 : a(9)=2008683 : a(10)=659226
a(11)=14709 : a(12)=5147 : a(13)=8859 : a(14)=2926 : a(15)=13
:
: REM Vector storage of equations' right-hand side.
:
x(1)=5937938 : x(2)=2046485 : x(3)=3526413 : x(4)=1130177 : x(5)=5003
:
FOR j% = 1 TO n% : REM Computing the decomposition of the normal equations.
q%=((j%*(j%+1)) DIV 2)
IF j%>1 THEN
FOR i% = j% TO n%
m% = ((i%*(i%-1)) DIV 2)+j% : s=a(m%)
FOR k% = 1 TO (j%-1) : s=s-a(m%-k%)*a(q%-k%) : NEXT k%
a(m%)=s
NEXT i%
ENDIF
: REM Cholesky factor.
s=SQR(a(q%))
FOR i% = j% TO n%
m%=((i%*(i%-1)) DIV 2)+j%
a(m%)=a(m%)/s
NEXT i%
PRINT
NEXT j%
:
VDU 30 : CLS : PRINT SPC(5);"Cholesky factor output vector, in row order."
FOR i% = 1 TO n%
FOR j% = 1 TO i%
PRINT SPC(5);a(j%)
NEXT j%
PRINT
NEXT i%
: REM Forward substitution commences.
x(1) = x(1)/a(1)
q%=1
FOR i%=2 TO n%
FOR j%=1 TO (i%-1)
q%=q%+1 : x(i%)=x(i%)-a(q%)*x(j%)
NEXT j%
q%=q%+1x(i%)=x(i%)/a(q%)
NEXT i%
: REM Backward substitution from here.
x(n%)=x(n%)/a(n%*(n%+1)/2)
FOR i%=n% TO 2 STEP -1
q%=i%*(i%-1)/2
FOR j%=1 TO (i%-1)
x(j%)=x(j%)-x(i%)*a(q%+j%)
NEXT j%
x(i%-1)=x(i%-1)/a(q%)
NEXT i%
:
PRINT SPC(5);"Determinate solution for the 5 unknowns in aproximately 30 milliseconds." : PRINT
FOR i% = 1 TO n%
PRINT SPC(5);"x(";STR$(i%);") = ";FNusing("###.##################",x(i%))
NEXT i%
:
:
END!!
:
:
|
|
 Logged
|
|
|
|
Richard Russell
Administrator
member is offline
Posts: 803
|
|
Re: Cholesky solution of matrices
« Reply #3 on: Jun 11th, 2017, 3:49pm » |
|
on Jun 11th, 2017, 1:56pm, Himself wrote:| REM The advent of the *FLOAT 80 function enhances the validity of the "most likely" solution |
|
Just to note that *FLOAT 80 has no effect at all in this program, since it affects only floating-point indirection and data files, and the program uses neither.
Richard.
|
|
Logged
|
|
|
|
|
