1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
! -*- f90 -*-
!
! Contains wrappers for the following LAPACK routines:
!
! Computational routines for the non-symmetric eigenproblem:
!
! gehrd (general Hessenberg reduction)
! gebal (general balancing)
! gebak (general backtransforming)
! orghr, unghr (orthogonal/unitary generate matrix after Hessenberg reduction) - Not Implemented
! ormhr, unmhr ((orthogonal/unitary multiply matrix after Hessenberg reduction) - Not Implemented
! hseqr (Hessenberg Schur factorization) - Not Implemented
! hsein (Hessenberg eigenvectors by inverse iteration) - Not Implemented
! trevc ((quasi)triangular eigenvectors) - Not Implemented
! trexc ((quasi)triangular reordering Schur factorization) - Not Implemented
! trsyl ((quasi)triangular Sylvester equation) - Not Implemented
! trsna ((quasi)triangular condition numbers of eigenvalues/vectors) - Not Implemented
! trsen ((quasi)triangular condition numbers of eigenvalue cluster/invariant subspace) - Not Implemented
!
subroutine <prefix>gehrd(n,lo,hi,a,tau,work,lwork,info)
!
! hq,tau,info = gehrd(a,lo=0,hi=n-1,lwork=n,overwrite_a=0)
! Reduce general matrix A to upper Hessenberg form H by unitary similarity
! transform Q^H * A * Q = H
!
! Q = H(lo) * H(lo+1) * ... * H(hi-1)
! H(i) = I - tau * v * v^H
! v[0:i+1] = 0, v[i+1]=1, v[hi+1:n] = 0
! v[i+2:hi+1] is stored in hq[i+2:hi+i,i]
! tau is tau[i]
!
! hq for n=7,lo=1,hi=5:
! [a a h h h h a
! a h h h h a
! h h h h h h
! v2h h h h h
! v2v3h h h h
! v2v3v4h h h
! a]
!
callstatement { hi++; lo++; (*f2py_func)(&n,&lo,&hi,a,&n,tau,work,&lwork,&info); }
callprotoargument int*,int*,int*,<ctype>*,int*,<ctype>*,<ctype>*,int*,int*
integer intent(hide),depend(a) :: n = shape(a,0)
<ftype> dimension(n,n),intent(in,out,copy,out=ht),check(shape(a,0)==shape(a,1)) :: a
integer intent(in),optional :: lo = 0
integer intent(in),optional,depend(n) :: hi = n-1
<ftype> dimension(n-1),intent(out),depend(n) :: tau
<ftype> dimension(lwork),intent(cahce,hide),depend(lwork) :: work
integer intent(in),optional,depend(n),check(lwork>=MAX(n,1)) :: lwork = MAX(n,1)
integer intent(out) :: info
end subroutine <prefix>gehrd
subroutine <prefix>gebal(scale,permute,n,a,m,lo,hi,pivscale,info)
!
! ba,lo,hi,pivscale,info = gebal(a,scale=0,permute=0,overwrite_a=0)
! Balance general matrix a.
! hi,lo are such that ba[i][j]==0 if i>j and j=0...lo-1 or i=hi+1..n-1
! pivscale([0:lo], [lo:hi+1], [hi:n+1]) = (p1,d,p2) where (p1,p2)[j] is
! the index of the row and column interchanged with row and column j.
! d[j] is the scaling factor applied to row and column j.
! The order in which the interchanges are made is n-1 to hi+1, then 0 to lo-1.
!
! P * A * P = [[T1,X,Y],[0,B,Z],[0,0,T2]]
! BA = [[T1,X*D,Y],[0,inv(D)*B*D,ind(D)*Z],[0,0,T2]]
! where D = diag(d), T1,T2 are upper triangular matrices.
! lo,hi mark the starting and ending columns of submatrix B.
callstatement { (*f2py_func)((permute?(scale?"B":"P"):(scale?"S":"N")),&n,a,&m,&lo,&hi,pivscale,&info); hi--; lo--; }
callprotoargument char*,int*,<ctype>*,int*,int*,int*,<ctypereal>*,int*
integer intent(in),optional :: permute = 0
integer intent(in),optional :: scale = 0
integer intent(hide),depend(a,n) :: m = shape(a,0)
integer intent(hide),depend(a) :: n = shape(a,1)
check(m>=n) m
integer intent(out) :: hi,lo
<ftypereal> dimension(n),intent(out),depend(n) :: pivscale
<ftype> dimension(m,n),intent(in,out,copy,out=ba) :: a
integer intent(out) :: info
end subroutine <prefix>gebal
|
