F
fogh
Guest
Hi All,
lately I used the functions below to test a few things about square
root of definite positive symmetric square matrices. Don t expect it
would know the socks of BLAS/LAPACK, but it did the job for me, even if
these are just naive SKILL implementation of numerical recipes. They
give/take the same types as the cadence ^matrix functions. The Jacobi
would be nicer if it returned a DPL.
;; MCreate
procedure( MCreate( n )
let((a V)
declare( a[n] )
for( i 0 n-1
a=declare(V[n])
for( j 0 n-1
a[j]=0.0
);j
);i
a ))
;; MIdent
procedure( MIdent
let((M)
M=MCreate
for( i 0 n - 1 M=1.0 )
M
))
;; MCopy
procedure( MCopy(a)
let((M V n)
n=length(a)
declare(M[n])
for( i 0 n-1
M=declare(V[n])
for( j 0 n-1
M[j]=a[j]
);j
);i
M
))
;; MSum
procedure( MSum(a b)
let((n M)
n=length(a)
M=MCreate
for( i 0 n-1
for( j 0 n-1
M[j]= a[j] + b[j]
);j
);i
M
))
;; MSubstract
procedure( MSubstract(a b)
let((n M)
n=length(a)
M=MCreate
for( i 0 n-1
for( j 0 n-1
M[j]=a[j] - b[j]
);j
);i
M
))
;; MProduct
procedure( MProduct(a b)
let((n M)
n=length(a)
M=MCreate
for( i 0 n-1
for( j 0 n-1
for( k 0 n-1
M[j]=M[j] + a[k] * b[k][j]
);k
);j
);i
M
))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; MPrint
procedure( MPrint(a)
let(
n=length(a)
for( i 0 n-1
for( j 0 n-1
printf("%10.04g" a[j] )
);j
printf("\n")
);i
));procedure
;; VPrint
procedure( VPrint(a)
let(
n=length(a)
for( i 0 n-1 printf("%10.04g" a )
printf("\n")
);i
));procedure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure( MTranspose(a)
let((n M)
n=length(a)
M=MCreate
for( i 0 n-1
for( j 0 n-1
M[j]=a[j]
);j
);i
M
));procedure
procedure( VDiag(a)
let((n v)
n=length(a)
declare(v[n])
for( i 0 n-1 v=a )
v
));procedure
procedure( MDiag(v)
let((n M)
n=length(v)
M=MCreate
for( i 0 n-1 M=v )
M
));procedure
;; MVector
procedure( MColVector(v)
let((M n)
n=length(v)
M=MCreate
for( i 0 n-1 M[0]=v )
M
))
;; MVProduct
procedure( MVProduct(a v)
let((n M)
n=length(a)
declare(M[n])
for( i 0 (n - 1)
M=0
for( j 0 (n - 1)
M=M+a[j]*v[j]
)
)
M
))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure( NRCcholdc(a)
/* Given a positive-definite symmetric matrix A, this routine constructs
its Cholesky decomposition, A = L ˇ LT .
The Cholesky factor L is returned */
let(
(i j k n;
sum;
l p
);local vars
n=length(a)
l=MCopy(a)
declare(p[n])
for(i 0 n-1
for(j i n-1
sum=l[j]
for( k 0 i-1 sum=sum-l[k]*l[j][k] )
if(i==j
then
unless(sum > 0.0 error("Choleski decomposition failed. Matrix is not
positive-definite ?") )
p=sqrt(sum)
else
l[j]=sum/p
);if i==j
);for j
);for i
for(i 0 n-1
l=p
for(j i+1 n-1
l[j]=0.0
)
)
l
);let
);procedure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure(NRCran1()
random(fix(1G))*random(fix(1G))/(1.0*1G**2)
)
procedure( NRCgasdev()
/* Returns a normally distributed deviate with zero mean and unit
variance, using ran1(idum) as the uniform deviates.*/
let(
( gset fac (rsq 0.0) v1 v2 )
while( rsq >= 1.0 || rsq == 0.0
;printf(".")
v1=2.0*NRCran1()-1.0; pick two uniform numbers in the square extending
from -1 to +1 in each direction,
v2=2.0*NRCran1()-1.0;
rsq=v1*v1+v2*v2 ; see if they are in the unit circle,
)
fac=sqrt(-2.0*log(rsq)/rsq); Now make the Box-Muller transformation to
get two normal deviates. Return one and save the other for next time.
gset=v1*fac;
))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure( VOne( n)
let((v)
declare( v[n] )
for( i 0 n-1
v=1.0
)
v
))
procedure(NRCrotate(a i j k l)
let( (g h)
g=a[j]
h=a[k][l]
a[j]=g-s*(h+g*tau)
a[k][l]=h+s*(g-h*tau)
))
procedure(NRCjacobi(A) ;;, int n, float d[], float **v, int *nrot
/*
Computes all eigenvalues and eigenvectors of a real symmetric matrix
a[1..n][1..n] . On
output, elements of a above the diagonal are destroyed. d[1..n] returns
the eigenvalues of a.
v[1..n][1..n] is a matrix whose columns contain, on output, the
normalized eigenvectors of
a. nrot returns the number of Jacobi rotations that were required.
*/
prog(
( j iq ip i;
tresh theta tau T sm s h g c
b z;
n d v nrot
maxiter
a
);local vars
a=MCopy(A)
maxiter=50
n=length(a)
d=VOne
v=MIdent ;;Initialize to the identity matrix.
b=VOne
z=VOne
;; Initialize b and d to the diagonal of a.
for(ip 0 n-1
b[ip]=d[ip]=a[ip][ip] ; This vector will accumulate terms
z[ip]=0.0 ; of the form tapq as in equation (11.1.14).
);for ip
nrot=0
for(i 0 maxiter-1
sm=0.0 ;Sum off-diagonal elements.
for(ip 0 n-2
for(iq ip+1 n-1 sm=sm+abs(a[ip][iq]) )
);for ip
if(sm == 0.0 ;The normal return
return(list(v d nrot))
);fi convergence
if(i<4
then ;...on the first three sweeps.
tresh=0.2*sm/(n*n)
else ;...thereafter.
tresh=0.0
);fi thresh
for(ip 0 n-2
for(iq ip+1 n-1
g=100.0*abs(a[ip][iq]) ;After four sweeps, skip the rotation if the
off-diagonal element is small.
if(i > 4 && (abs(d[ip])+g) == abs(d[ip])
&& (abs(d[iq])+g) == abs(d[iq])
then a[ip][iq]=0.0
else if(abs(a[ip][iq]) > tresh
then
h=d[iq]-d[ip]
if((abs(h)+g) == abs(h) ;t = 1/(2 theta)
then T=a[ip][iq]/h
else ; Equation (11.1.10).
theta=0.5*h/(a[ip][iq])
T=1.0/(abs(theta)+sqrt(1.0+theta*theta))
if(theta < 0.0 then T=-T )
);fi abs(h)+g
c=1.0/sqrt(1+T*T)
s=T*c
tau=s/(1.0+c)
h=T*a[ip][iq]
z[ip]=z[ip]-h
z[iq]=z[iq]+h
d[ip]=d[ip]-h
d[iq]=d[iq]+h
a[ip][iq]=0.0 ;Case of rotations 1 j < p.
for(j 0 ip-1 NRCrotate(a j ip j iq));Case of rotations p < j < q .
for(j ip+1 iq-1 NRCrotate(a ip j j iq));Case of rotations q < j n.
for(j iq+1 n-1 NRCrotate(a ip j iq j))
for(j 0 n-1 NRCrotate(v j ip j iq))
++nrot
);fi abs(a[ip][iq]) > tresh
);fi i>4 ...
);for iq
);for ip
for(ip 0 n-1
b[ip]=b[ip]+ z[ip] ;Update d with the sum of tapq ,
d[ip]=b[ip] ;and reinitialize z.
z[ip]=0.0
);for ip
); for i to maxiter
error("Too many iterations in routine jacobi");
);let
);procedure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure(MRandomSymetric
let((a V)
declare( a[n] )
for( i 0 n-1 a=declare(V[n]) )
for( i 0 n-1
for( j i n-1
a[j]=a[j]=NRCran1()
);j
);i
a ))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure(Vfunc(v func)
let((r n)
n=length(v)
declare( r[n] )
for( i 0 n-1
r=funcall(func v)
)
r
))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure( Mradius( a )
let((n r)
n=length(a)
r=0.0
for( i 0 n-1
for( j 0 n-1
r=r+abs(a[j])**2
);j
);i
sqrt(r)
))
;c=MRandomSymetric(size)
el=NRCjacobi(c)
printf("eigenvectors:\n") MPrint(car(el))
printf("eigenvalues: \n") VPrint(cadr(el))
printf("Jacobi rotations: \n") printf("%L\n" caddr(el))
V=car(el) ;eigenvectors
L=cadr(el) ;eigenvalues
ML=MDiag(L)
VsqrL=Vfunc(L 'sqrt) ;square root of eigenvalues ;;rem: small neg values
gives complex.
MsqrL=MDiag(VsqrL) ;diag matrix of square root of eigenvalues
/* check Jacobi precision.
println(Mradius(
MSubstract(
MProduct(V MProduct(ML MTranspose(V)))
c
)
))
;=> 1E-16
println(Mradius(
MSubstract(
MProduct(
MProduct(V MsqrL)
MTranspose(MProduct(V MsqrL))
);prod
c
);subs
));print
;=> 1E-16
*/
;MPrint(MProduct(V MProduct(MsqrL MTranspose(V)) ) )
a=MProduct(V MProduct(MsqrL MTranspose(V)) )
println(Mradius(
MSubstract(
MProduct(MTranspose(a) a)
c
);-
));print
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; a is symetric , and radius(c-aT.a) <10**-15
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
lately I used the functions below to test a few things about square
root of definite positive symmetric square matrices. Don t expect it
would know the socks of BLAS/LAPACK, but it did the job for me, even if
these are just naive SKILL implementation of numerical recipes. They
give/take the same types as the cadence ^matrix functions. The Jacobi
would be nicer if it returned a DPL.
;; MCreate
procedure( MCreate( n )
let((a V)
declare( a[n] )
for( i 0 n-1
a=declare(V[n])
for( j 0 n-1
a[j]=0.0
);j
);i
a ))
;; MIdent
procedure( MIdent
let((M)
M=MCreate
for( i 0 n - 1 M=1.0 )
M
))
;; MCopy
procedure( MCopy(a)
let((M V n)
n=length(a)
declare(M[n])
for( i 0 n-1
M=declare(V[n])
for( j 0 n-1
M[j]=a[j]
);j
);i
M
))
;; MSum
procedure( MSum(a b)
let((n M)
n=length(a)
M=MCreate
for( i 0 n-1
for( j 0 n-1
M[j]= a[j] + b[j]
);j
);i
M
))
;; MSubstract
procedure( MSubstract(a b)
let((n M)
n=length(a)
M=MCreate
for( i 0 n-1
for( j 0 n-1
M[j]=a[j] - b[j]
);j
);i
M
))
;; MProduct
procedure( MProduct(a b)
let((n M)
n=length(a)
M=MCreate
for( i 0 n-1
for( j 0 n-1
for( k 0 n-1
M[j]=M[j] + a[k] * b[k][j]
);k
);j
);i
M
))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; MPrint
procedure( MPrint(a)
let(
n=length(a)
for( i 0 n-1
for( j 0 n-1
printf("%10.04g" a[j] )
);j
printf("\n")
);i
));procedure
;; VPrint
procedure( VPrint(a)
let(
n=length(a)
for( i 0 n-1 printf("%10.04g" a )
printf("\n")
);i
));procedure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure( MTranspose(a)
let((n M)
n=length(a)
M=MCreate
for( i 0 n-1
for( j 0 n-1
M[j]=a[j]
);j
);i
M
));procedure
procedure( VDiag(a)
let((n v)
n=length(a)
declare(v[n])
for( i 0 n-1 v=a )
v
));procedure
procedure( MDiag(v)
let((n M)
n=length(v)
M=MCreate
for( i 0 n-1 M=v )
M
));procedure
;; MVector
procedure( MColVector(v)
let((M n)
n=length(v)
M=MCreate
for( i 0 n-1 M[0]=v )
M
))
;; MVProduct
procedure( MVProduct(a v)
let((n M)
n=length(a)
declare(M[n])
for( i 0 (n - 1)
M=0
for( j 0 (n - 1)
M=M+a[j]*v[j]
)
)
M
))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure( NRCcholdc(a)
/* Given a positive-definite symmetric matrix A, this routine constructs
its Cholesky decomposition, A = L ˇ LT .
The Cholesky factor L is returned */
let(
(i j k n;
sum;
l p
);local vars
n=length(a)
l=MCopy(a)
declare(p[n])
for(i 0 n-1
for(j i n-1
sum=l[j]
for( k 0 i-1 sum=sum-l[k]*l[j][k] )
if(i==j
then
unless(sum > 0.0 error("Choleski decomposition failed. Matrix is not
positive-definite ?") )
p=sqrt(sum)
else
l[j]=sum/p
);if i==j
);for j
);for i
for(i 0 n-1
l=p
for(j i+1 n-1
l[j]=0.0
)
)
l
);let
);procedure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure(NRCran1()
random(fix(1G))*random(fix(1G))/(1.0*1G**2)
)
procedure( NRCgasdev()
/* Returns a normally distributed deviate with zero mean and unit
variance, using ran1(idum) as the uniform deviates.*/
let(
( gset fac (rsq 0.0) v1 v2 )
while( rsq >= 1.0 || rsq == 0.0
;printf(".")
v1=2.0*NRCran1()-1.0; pick two uniform numbers in the square extending
from -1 to +1 in each direction,
v2=2.0*NRCran1()-1.0;
rsq=v1*v1+v2*v2 ; see if they are in the unit circle,
)
fac=sqrt(-2.0*log(rsq)/rsq); Now make the Box-Muller transformation to
get two normal deviates. Return one and save the other for next time.
gset=v1*fac;
))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure( VOne( n)
let((v)
declare( v[n] )
for( i 0 n-1
v=1.0
)
v
))
procedure(NRCrotate(a i j k l)
let( (g h)
g=a[j]
h=a[k][l]
a[j]=g-s*(h+g*tau)
a[k][l]=h+s*(g-h*tau)
))
procedure(NRCjacobi(A) ;;, int n, float d[], float **v, int *nrot
/*
Computes all eigenvalues and eigenvectors of a real symmetric matrix
a[1..n][1..n] . On
output, elements of a above the diagonal are destroyed. d[1..n] returns
the eigenvalues of a.
v[1..n][1..n] is a matrix whose columns contain, on output, the
normalized eigenvectors of
a. nrot returns the number of Jacobi rotations that were required.
*/
prog(
( j iq ip i;
tresh theta tau T sm s h g c
b z;
n d v nrot
maxiter
a
);local vars
a=MCopy(A)
maxiter=50
n=length(a)
d=VOne
v=MIdent
;;Initialize to the identity matrix.b=VOne
z=VOne
;; Initialize b and d to the diagonal of a.
for(ip 0 n-1
b[ip]=d[ip]=a[ip][ip] ; This vector will accumulate terms
z[ip]=0.0 ; of the form tapq as in equation (11.1.14).
);for ip
nrot=0
for(i 0 maxiter-1
sm=0.0 ;Sum off-diagonal elements.
for(ip 0 n-2
for(iq ip+1 n-1 sm=sm+abs(a[ip][iq]) )
);for ip
if(sm == 0.0 ;The normal return
return(list(v d nrot))
);fi convergence
if(i<4
then ;...on the first three sweeps.
tresh=0.2*sm/(n*n)
else ;...thereafter.
tresh=0.0
);fi thresh
for(ip 0 n-2
for(iq ip+1 n-1
g=100.0*abs(a[ip][iq]) ;After four sweeps, skip the rotation if the
off-diagonal element is small.
if(i > 4 && (abs(d[ip])+g) == abs(d[ip])
&& (abs(d[iq])+g) == abs(d[iq])
then a[ip][iq]=0.0
else if(abs(a[ip][iq]) > tresh
then
h=d[iq]-d[ip]
if((abs(h)+g) == abs(h) ;t = 1/(2 theta)
then T=a[ip][iq]/h
else ; Equation (11.1.10).
theta=0.5*h/(a[ip][iq])
T=1.0/(abs(theta)+sqrt(1.0+theta*theta))
if(theta < 0.0 then T=-T )
);fi abs(h)+g
c=1.0/sqrt(1+T*T)
s=T*c
tau=s/(1.0+c)
h=T*a[ip][iq]
z[ip]=z[ip]-h
z[iq]=z[iq]+h
d[ip]=d[ip]-h
d[iq]=d[iq]+h
a[ip][iq]=0.0 ;Case of rotations 1 j < p.
for(j 0 ip-1 NRCrotate(a j ip j iq));Case of rotations p < j < q .
for(j ip+1 iq-1 NRCrotate(a ip j j iq));Case of rotations q < j n.
for(j iq+1 n-1 NRCrotate(a ip j iq j))
for(j 0 n-1 NRCrotate(v j ip j iq))
++nrot
);fi abs(a[ip][iq]) > tresh
);fi i>4 ...
);for iq
);for ip
for(ip 0 n-1
b[ip]=b[ip]+ z[ip] ;Update d with the sum of tapq ,
d[ip]=b[ip] ;and reinitialize z.
z[ip]=0.0
);for ip
); for i to maxiter
error("Too many iterations in routine jacobi");
);let
);procedure
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure(MRandomSymetric
let((a V)
declare( a[n] )
for( i 0 n-1 a=declare(V[n]) )
for( i 0 n-1
for( j i n-1
a[j]=a[j]=NRCran1()
);j
);i
a ))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure(Vfunc(v func)
let((r n)
n=length(v)
declare( r[n] )
for( i 0 n-1
r=funcall(func v)
)
r
))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
procedure( Mradius( a )
let((n r)
n=length(a)
r=0.0
for( i 0 n-1
for( j 0 n-1
r=r+abs(a[j])**2
);j
);i
sqrt(r)
))
;c=MRandomSymetric(size)
el=NRCjacobi(c)
printf("eigenvectors:\n") MPrint(car(el))
printf("eigenvalues: \n") VPrint(cadr(el))
printf("Jacobi rotations: \n") printf("%L\n" caddr(el))
V=car(el) ;eigenvectors
L=cadr(el) ;eigenvalues
ML=MDiag(L)
VsqrL=Vfunc(L 'sqrt) ;square root of eigenvalues ;;rem: small neg values
gives complex.
MsqrL=MDiag(VsqrL) ;diag matrix of square root of eigenvalues
/* check Jacobi precision.
println(Mradius(
MSubstract(
MProduct(V MProduct(ML MTranspose(V)))
c
)
))
;=> 1E-16
println(Mradius(
MSubstract(
MProduct(
MProduct(V MsqrL)
MTranspose(MProduct(V MsqrL))
);prod
c
);subs
;=> 1E-16
*/
;MPrint(MProduct(V MProduct(MsqrL MTranspose(V)) ) )
a=MProduct(V MProduct(MsqrL MTranspose(V)) )
println(Mradius(
MSubstract(
MProduct(MTranspose(a) a)
c
);-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; a is symetric , and radius(c-aT.a) <10**-15
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;