! 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