! examples of matrix operations
! =============================
USE cbmat
USE rtchecks

COMPLEX M(0..4,0..4),t(0..4),x(0..4)
DO M(i,j)=(1.+j*I)**(i) FOR ALL i,j
t=0; t(4)=24 ! initialize some arrays and matrices

TYPEOF(M(3,*)) temp=M(3,*)
M(3,*)=M(0,*)
M(0,*)=temp ! swap 0th and 3rd row

Lanczos(x,M,t,1E-10) ! solve linear system by Lanczos method
WRITE x
WRITE ABS(t-M*x) ! and test the result

plu=PLU(M) ! Perform LU decomposition with pivoting
x=plu\t    ! and apply it to solve a linear system
WRITE x    ! This is faster than storing the inverse.
WRITE ABS(t-M*x) ! Test the result

TYPEOF(M) M1,M2  ! prepare for copies
M1=M;LUdecomp M1 ! LU decomposition (no pivoting) in place
x=M1\t ! division to the left, with no stored inverse
WRITE ABS(t-M*x) ! and test
x=t/M1 ! division to the right
WRITE ABS(t-x*M) ! and test
DO M2(i,*)=M(i,*)/plu FOR ALL i; WRITE M2 ! divide each row, equivalent to
                                   ! right multiplication with the inverse

M1=INV(M) ! M1 is the inverse matrix

M2($i,$j)=CONJG(M1($i,$k))*M($k,$j); WRITE M2
        ! matrix product in Einstein notation

STORED ARRAY(0..4,-2..2) OF COMPLEX Ms ! Same operations can be performed
                                       ! on STORED arrays
Ms(*,*)=M(*,*+2);LUdecomp Ms ! LU decomposition on disk
x=Ms\t                       ! right division
WRITE ABS(t-M(*,*+2)*x)      ! test
x=t/Ms                       ! left division
WRITE ABS(t-x*M(*,*+2))      ! test
Ms(*,*)=M(*,*+2);LUdecomp M(*,*+2)
x=M(*,*+2)\t      ! opposite roles
WRITE ABS(t-Ms*x)
x=t/M(*,*+2)
WRITE ABS(t-x*Ms)

ARRAY(1..5,1..5) OF COMPLEX A=1,B=1,C

C=A*B ! direct matrix multiplication (use sparingly!)
LUdecomp A
B=A\C ! and division