jashwanth9/parrot-lapack forked from leto/parrot-lapack

Documentation added

• Loading branch information...
1 parent 5389325 commit 0baee2ca63246a33bd89b5afbd2cc5edeac8c73f committed Sep 5, 2012
56 src/LAPACK/C_Double.winxed
 @@ -1,7 +1,14 @@ +// ZGETRF computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges. + namespace zgetrf_func{ const int PRINT_DEBUG_STUFF = 0; +/* + A: COMPLEX*16 array, dimension (LDA,N) + IPIV: INTEGER array, dimension (min(M,N)) +*/ + function zgetrf_exec(var A,int ipiv) { var pla = loadlib("linalg_group"); @@ -12,6 +19,14 @@ function zgetrf_exec(var A,int ipiv) say("Given Matrix:"); say(A); +/* + A is COMPLEX*16 array, dimension (LDA,N) + On entry, the M-by-N matrix to be factored. + On exit, the factors L and U from the factorization + A = P*L*U; the unit diagonal elements of L are not stored. + */ + + int m,n,lda,ipiv_size,info; m=A.rows; n=A.cols; @@ -25,6 +40,16 @@ function zgetrf_exec(var A,int ipiv) zgetrf(m,n,A,lda,ipiv,info); return info; +/* +INFO is INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value + > 0: if INFO = i, U(i,i) is exactly zero. The factorization + has been completed, but the factor U is exactly + singular, and division by zero will occur if it is used + to solve a system of equations. +*/ + } function max(var a,var b) @@ -56,34 +81,3 @@ function debug(var matrix, string msg, var args [slurpy]) } } - -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'ComplexMatrix2D'; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - using zgetrf_func.ipiv_size; - int ipiv_size_v; - ipiv_size_v=ipiv_size(a); - say(ipiv_size_v); - int ipiv[ipiv_size_v]; - int i; - for(i=0;i
62 src/LAPACK/C_eigenvalues.winxed
 @@ -1,13 +1,33 @@ +//ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors. +//The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. + + namespace zgeev_func_eva{ const int PRINT_DEBUG_STUFF = 0; +/* + A: COMPLEX*16 array, dimension (LDA,N) + W: COMPLEX*16 array, dimension (N) + JOBVL: CHARACTER + JOBVR: CHARACTER + = 'N': left eigenvectors of A are not computed; + = 'V': left eigenvectors of A are computed. +*/ + function zgeev_exec_eva(var A,var W,int jobvl,int jobvr) { var pla = loadlib("linalg_group"); var lapack = loadlib('liblapack.so'); if (lapack == null || !lapack) die("Cannot find liblapack. Search for the correct paths"); + +/* + A is COMPLEX*16 array, dimension (LDA,N) + On entry, the N-by-N matrix A. + On exit, A has been overwritten. +*/ + var VL=new 'ComplexMatrix2D'; var VR=new 'ComplexMatrix2D'; int m,n,lda,ldvl,ldvr,lwork,info,rwork; @@ -35,7 +55,19 @@ function zgeev_exec_eva(var A,var W,int jobvl,int jobvr) zgeev(jobvl,jobvr,n,A,lda,W,VL,ldvl,VR,ldvr,WORK,lwork,rwork,info); + return info; + +/* + INFO is INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value. + > 0: if INFO = i, the QR algorithm failed to compute all the + eigenvalues, and no eigenvectors have been computed; + elements i+1:N of WR and WI contain eigenvalues which + have converged. + */ + } function max(var a,var b) @@ -44,33 +76,3 @@ function max(var a,var b) } } - -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'ComplexMatrix2D'; - var w=new 'ComplexMatrix2D'; - - int jobvl,jobvr; - jobvl=jobvr= ord('N', 0); - var info; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - - using zgeev_func_eva.zgeev_exec_eva; - info=zgeev_exec_eva(a,w,jobvl,jobvr); - - if(info==0) - { - say("successful\n"); - say("eigen values=\n"); - say(w); - } -} - -*/
74 src/LAPACK/C_eigenvectors.winxed
 @@ -1,7 +1,20 @@ +//ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors. +//The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real. + + namespace zgeev_func_eve{ const int PRINT_DEBUG_STUFF = 0; +/* + A: COMPLEX*16 array, dimension (LDA,N) + W: COMPLEX*16 array, dimension (N) + JOBVL: CHARACTER + JOBVR: CHARACTER + = 'N': left eigenvectors of A are not computed; + = 'V': left eigenvectors of A are computed. +*/ + function zgeev_exec_eve(var A,var W,var VL,var VR) { var pla = loadlib("linalg_group"); @@ -16,6 +29,12 @@ function zgeev_exec_eve(var A,var W,var VL,var VR) if(zgeev == null || !zgeev) die("Not zgeev"); +/* + A is COMPLEX*16 array, dimension (LDA,N) + On entry, the N-by-N matrix A. + On exit, A has been overwritten. +*/ + var WORK=new 'ComplexMatrix2D'; jobvl=jobvr= ord('V', 0); m=A.rows; @@ -33,9 +52,36 @@ function zgeev_exec_eve(var A,var W,var VL,var VR) WORK.resize(lwork,lwork); } +/* +VL is COMPLEX*16 array, dimension (LDVL,N) + If JOBVL = 'V', the left eigenvectors u(j) are stored one + after another in the columns of VL, in the same order + as their eigenvalues. + If JOBVL = 'N', VL is not referenced. + u(j) = VL(:,j), the j-th column of VL. + +VR is COMPLEX*16 array, dimension (LDVR,N) + If JOBVR = 'V', the right eigenvectors v(j) are stored one + after another in the columns of VR, in the same order + as their eigenvalues. + If JOBVR = 'N', VR is not referenced. + v(j) = VR(:,j), the j-th column of VR. + +*/ zgeev(jobvl,jobvr,n,A,lda,W,VL,ldvl,VR,ldvr,WORK,lwork,rwork,info); + return info; +/* + INFO is INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value. + > 0: if INFO = i, the QR algorithm failed to compute all the + eigenvalues, and no eigenvectors have been computed; + elements i+1:N of WR and WI contain eigenvalues which + have converged. + */ + } function max(var a,var b) @@ -44,31 +90,3 @@ function max(var a,var b) } } - -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'ComplexMatrix2D'; - var w=new 'ComplexMatrix2D'; - var vl=new 'ComplexMatrix2D'; - var vr=new 'ComplexMatrix2D'; - var info; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - - using zgeev_func_eve.zgeev_exec_eve; - info=zgeev_exec_eve(a,w,vl,vr); - - if(info==0) - { - say("right eigenvectors=\n",vr); - say("left eigenvectors=\n",vl); - } -} - -*/
44 src/LAPACK/C_inverse.winxed
 @@ -1,9 +1,15 @@ +//ZGETRI computes the inverse of a matrix using the LU factorization computed by ZGETRF. +//This method inverts U and then computes inv(A) by solving the system inv(A)*L = inv(U) for inv(A). + + \$load "C_Double.pbc"; namespace zgetri_func{ const int PRINT_DEBUG_STUFF = 0; +// A is COMPLEX*16 array, dimension (LDA,N) + function zgetri_exec(var a) { var pla = loadlib("linalg_group"); @@ -21,8 +27,15 @@ function zgetri_exec(var a) ipiv[i]=0; say("Given Matrix:"); - say(a); + say(a); +/* + A is COMPLEX*16 array, dimension (LDA,N) + On entry, the factors L and U from the factorization + A = P*L*U as computed by ZGETRF. + On exit, if INFO = 0, the inverse of the original matrix A. +*/ + var zgetri = dlfunc(lapack, "zgetri_", "vppppppp"); if(zgetri == null || !zgetri) die("Not ZGETRI"); @@ -42,6 +55,15 @@ function zgetri_exec(var a) zgetri(n,a,lda,ipiv,work,lwork,info); return info; +/* +INFO is INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value + > 0: if INFO = i, U(i,i) is exactly zero; the matrix is + singular and its inverse could not be computed. + +*/ + } @@ -66,23 +88,3 @@ function debug(var matrix, string msg, var args [slurpy]) } -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'ComplexMatrix2D'; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - using zgetri_func.zgetri_exec; - int info; - info=zgetri_exec(a); - if(info==0) - say(a); - say(""); -} - -*/
62 src/LAPACK/C_linearequation.winxed
 @@ -1,7 +1,27 @@ +/* +ZSPSV computes the solution to a complex system of linear equations + A * X = B, + where A is an N-by-N symmetric matrix stored in packed format and X and B are N-by-NRHS matrices. + + The diagonal pivoting method is used to factor A as + A = U * D * U**T, if UPLO = 'U', or + A = L * D * L**T, if UPLO = 'L', + where U (or L) is a product of permutation and unit upper (lower) triangular matrices, D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then used to solve the system of equations A * X = B. + +*/ + namespace zspsv_func{ const int PRINT_DEBUG_STUFF = 0; +/* + A: DOUBLE PRECISION array, dimension (LDA,N) + B: DOUBLE PRECISION array, dimension (LDB,NRHS) + UPLO: CHARACTER + 'U': Upper triangle of A is stored + 'L': Lower triangle of A is stored. +*/ + function zspsv_exec(var A,var b,int uplo) { @@ -10,6 +30,20 @@ function zspsv_exec(var A,var b,int uplo) if (lapack == null || !lapack) die("Cannot find liblapack. Search for the correct paths"); +/* + AP is COMPLEX*16 array, dimension (N*(N+1)/2) + On entry, the upper or lower triangle of the symmetric matrix + A, packed columnwise in a linear array. The j-th column of A + is stored in the array AP as follows: + if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; + if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. + + On exit, the block diagonal matrix D and the multipliers used + to obtain the factor U or L from the factorization + A = U*D*U**T or A = L*D*L**T as computed by ZSPTRF, stored as + a packed triangular matrix in the same storage format as A. +*/ + int m,n,nrhs,ldb,info; m=A.rows; n=A.cols; @@ -84,31 +118,3 @@ function debug(var matrix, string msg, var args [slurpy]) } -/* - -function main[main](var args) -{ - var pla = loadlib("linalg_group"); - var a = new 'ComplexMatrix2D'; - var b=new 'ComplexMatrix2D'; - var ap=new 'ComplexMatrix2D'; - int uplo; - uplo=85; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - b.initialize_from_args(3, 1, - 1.0, - 1.0, - 1.0); - using zspsv_func.zspsv_exec; - ap=zspsv_exec(a,b,uplo); - - say("AP="); - say(ap); - say("B="); - say(b); -} - -*/
43 src/LAPACK/C_qrfactorisation.winxed
 @@ -1,7 +1,14 @@ +// ZGEQRF computes a QR factorization of a real M-by-N matrix A: +//A = Q * R. + namespace zgeqrf_func{ const int PRINT_DEBUG_STUFF = 0; +/* +A : COMPLEX*16 array, dimension (LDA,N) +*/ + function zgeqrf_exec(var A) { @@ -13,6 +20,17 @@ function zgeqrf_exec(var A) say("Given Matrix:"); say(A); +/* +A is COMPLEX*16 array, dimension (LDA,N) + On entry, the M-by-N matrix A. + On exit, the elements on and above the diagonal of the array + contain the min(M,N)-by-N upper trapezoidal matrix R (R is + upper triangular if m >= n); the elements below the diagonal, + with the array TAU, represent the orthogonal matrix Q as a + product of min(m,n) elementary reflectors (see Further + Details). +*/ + int m,n,lda,lwork,tau_s,work_s,info; m=A.rows; @@ -37,6 +55,11 @@ function zgeqrf_exec(var A) return info; +/* + INFO is INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value +*/ } @@ -60,23 +83,3 @@ function debug(var matrix, string msg, var args [slurpy]) } -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'ComplexMatrix2D'; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - using zgeqrf_func.zgeqrf_exec; - int info; - info=zgeqrf_exec(a); - if(info==0) - say(a); -} - - -*/
85 src/LAPACK/C_singlevaluedecomposition.winxed
 @@ -1,7 +1,33 @@ +/* +ZGESVD computes the singular value decomposition (SVD) of a complex M-by-N matrix A, optionally computing the left and/or right singular vectors. The SVD is written + A = U * SIGMA * conjugate-transpose(V) + + where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M unitary matrix, and V is an N-by-N unitary matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A. + + Note that the routine returns V**H, not V. +*/ + namespace zgesvd_func{ const int PRINT_DEBUG_STUFF = 0; +/* +A : COMPLEX*16 array, dimension (LDA,N) +JOBU : CHARACTER +JOBVT: CHARACTER + Specifies options for computing all or part of the matrix U: + = 'A': all M columns of U are returned in array U: + = 'S': the first min(m,n) columns of U (the left singular + vectors) are returned in the array U; + = 'O': the first min(m,n) columns of U (the left singular + vectors) are overwritten on the array A; + = 'N': no columns of U (no left singular vectors) are + computed. + S : DOUBLE PRECISION array, dimension (min(M,N)) + U : COMPLEX*16 array, dimension (LDU,UCOL) + VT : COMPLEX*16 array, dimension (LDVT,N) +*/ + function zgesvd_exec(var A,int jobu,int jobvt,var S,var U,var VT) { @@ -14,6 +40,20 @@ function zgesvd_exec(var A,int jobu,int jobvt,var S,var U,var VT) say("Given Matrix:"); say(A); +/* + A is COMPLEX*16 array, dimension (LDA,N) + On entry, the M-by-N matrix A. + On exit, + if JOBU = 'O', A is overwritten with the first min(m,n) + columns of U (the left singular vectors, + stored columnwise); + if JOBVT = 'O', A is overwritten with the first min(m,n) + rows of V**H (the right singular vectors, + stored rowwise); + if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A + are destroyed. +*/ + int m,n,lda,ldu,ldvt,lwork,info; m=A.rows; n=A.cols; @@ -46,6 +86,16 @@ function zgesvd_exec(var A,int jobu,int jobvt,var S,var U,var VT) zgesvd(jobu,jobvt,m,n,A,lda,S,U,ldu,VT,ldvt,work,lwork,rwork,info); return info; +/* +INFO is INTEGER + = 0: successful exit. + < 0: if INFO = -i, the i-th argument had an illegal value. + > 0: if ZBDSQR did not converge, INFO specifies how many + superdiagonals of an intermediate bidiagonal form B + did not converge to zero. See the description of RWORK + above for details. +*/ + } function max(var a,var b) @@ -68,38 +118,3 @@ function debug(var matrix, string msg, var args [slurpy]) } -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'ComplexMatrix2D'; - var S = new 'ComplexMatrix2D'; - var U = new 'ComplexMatrix2D'; - var VT = new 'ComplexMatrix2D'; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - - int jobu,jobvt; - jobu=jobvt= ord('A', 0); - - using zgesvd_func.zgesvd_exec; - int info; - - info=zgesvd_exec(a,jobu,jobvt,S,U,VT); - - if(info == 0) - { - say("successful"); - say(a); - say(S); - say(U); - say(VT); - } - -} - -*/
42 src/LAPACK/Double.winxed
 @@ -5,6 +5,11 @@ namespace dgetrf_func{ const int PRINT_DEBUG_STUFF = 0; +/* + A: DOUBLE PRECISION array, dimension (LDA,N) + IPIV: INTEGER array, dimension (min(M,N)) +*/ + function dgetrf_exec(var A,int ipiv) { var pla = loadlib("linalg_group"); @@ -16,6 +21,13 @@ function dgetrf_exec(var A,int ipiv) say("Given Matrix:"); say(A); + /* + A is DOUBLE PRECISION array, dimension (LDA,N) + On entry, the M-by-N matrix to be factored. + On exit, the factors L and U from the factorization + A = P*L*U; the unit diagonal elements of L are not stored. + */ + int m,n,lda,ipiv_size,info; m=A.rows; n=A.cols; @@ -105,33 +117,3 @@ function debug(var matrix, string msg, var args [slurpy]) } -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'NumMatrix2D'; - a.initialize_from_args(3, 3, - 2.0, -1.0, 1.0, - 4.0, 1.0, -1.0, - 1.0, 1.0, 1.0,); - using dgetrf_func.ipiv_size; - int ipiv_size_v; - ipiv_size_v=ipiv_size(a); - say(ipiv_size_v); - int ipiv[ipiv_size_v]; - int i; - for(i=0;i
27 src/LAPACK/LUdecomposition.winxed
 @@ -1,3 +1,5 @@ +//LUDECOMP_FUNC compputes the uform and Lform of a square matrix + \$load "inverse.pbc"; \$load "multiplication.pbc"; @@ -57,8 +59,6 @@ function Lform(var A,var u) using dgetri_func.dgetri_exec; int info; info=dgetri_exec(u); -//say("Here is the error:"); -// l=A*u; using multiplication.matmul; l=matmul(A,u); return l; @@ -71,21 +71,10 @@ function abs(var a) return (a); } -/* -function multiply(var P,var Q) -{ - - - - -} - -*/ - } - +/* function main[main](var args) { @@ -100,21 +89,15 @@ function main[main](var args) 1.0, 1.0, 1.0); using ludecompo_func.Uform; u=Uform(a); - say("Successful"); say("U="); say(u); say("A="); say(a); using ludecompo_func.Lform; l=Lform(a,u); - say("Successful"); say("L="); say(l); - -using dgetri_func.dgetri_exec; - int info; - info=dgetri_exec(u); - - } + +*/
79 src/LAPACK/eigenvalues.winxed
 @@ -1,16 +1,39 @@ +// DGEEV computes for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors. +// The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real + + namespace dgeev_func_eva{ const int PRINT_DEBUG_STUFF = 0; +/* + A: DOUBLE PRECISION array, dimension (LDA,N) + WR: DOUBLE PRECISION array, dimension (N) + WI: DOUBLE PRECISION array, dimension (N) + JOBVL: CHARACTER + JOBVR: CHARACTER + = 'N': left eigenvectors of A are not computed; + = 'V': left eigenvectors of A are computed. +*/ + function dgeev_exec_eva(var A,var WR,var WI,int jobvl,int jobvr) { var pla = loadlib("linalg_group"); var lapack = loadlib('liblapack.so'); if (lapack == null || !lapack) die("Cannot find liblapack. Search for the correct paths"); + + say("Given Matrix:"); + say(A); +/* +A is DOUBLE PRECISION array, dimension (LDA,N) + On entry, the N-by-N matrix A. + On exit, A has been overwritten. +*/ + var VL=new 'NumMatrix2D'; var VR=new 'NumMatrix2D'; - int m,n,lda,ldvl,ldvr,lwork,info; + int m,n,lda,ldvl,ldvr,lwork,info; var dgeev = dlfunc(lapack, "dgeev_", "vpppppppppppppp"); if(dgeev == null || !dgeev) @@ -21,8 +44,6 @@ function dgeev_exec_eva(var A,var WR,var WI,int jobvl,int jobvr) m=A.rows; n=A.cols; - if(m==n) - { lda=max(1,n); WR.resize(n,n); WI.resize(n,n); @@ -31,53 +52,29 @@ function dgeev_exec_eva(var A,var WR,var WI,int jobvl,int jobvr) VR.resize(ldvr,n); lwork=max(1,n*4); WORK.resize(lwork,lwork); - } dgeev(jobvl,jobvr,n,A,lda,WR,WI,VL,ldvl,VR,ldvr,WORK,lwork,info); + +/* + INFO is INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value. + > 0: if INFO = i, the QR algorithm failed to compute all the + eigenvalues, and no eigenvectors have been computed; + elements i+1:N of WR and WI contain eigenvalues which + have converged. + */ + return info; + } +//computes the max value + function max(var a,var b) { return a>b?a:b; } } - -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'NumMatrix2D'; - var wr=new 'NumMatrix2D'; - var wi=new 'NumMatrix2D'; - - int jobvl,jobvr; - jobvl=jobvr= ord('N', 0); - var info; - a.initialize_from_args(3, 3, - 11.0, 2.0, 3.0, - 4.0, 5.0, 6.0, - 7.0, 1.0, 2.0); - - using dgeev_func_eva.dgeev_exec_eva; - info=dgeev_exec_eva(a,wr,wi,jobvl,jobvr); - - if(info==0) - { - say("successful\n"); - say("eigen values=\n"); - int i,j; - for(i=0;i<3;++i) - { - for(j=0;j<3;++j) - print(wr[i,j],"+j",wi[i,j],"\t"); - say(); - } - } -} - -*/
77 src/LAPACK/eigenvectors.winxed
 @@ -1,75 +1,80 @@ +// DGEEV computes for an N-by-N real nonsymmetric matrix A, the eigenvalues and, optionally, the left and/or right eigenvectors. +// The computed eigenvectors are normalized to have Euclidean norm equal to 1 and largest component real +// This library is spefically fo calculating eigenvectors the left and/or right . + namespace dgeev_func_eve{ const int PRINT_DEBUG_STUFF = 0; +/* + A: DOUBLE PRECISION array, dimension (LDA,N) + WR: DOUBLE PRECISION array, dimension (N) + WI: DOUBLE PRECISION array, dimension (N) + VL: DOUBLE PRECISION array, dimension (LDVL,N) + VR: DOUBLE PRECISION array, dimension (LDVR,N) + +*/ + function dgeev_exec_eve(var A,var WR,var WI,var VL,var VR) { var pla = loadlib("linalg_group"); var lapack = loadlib('liblapack.so'); if (lapack == null || !lapack) die("Cannot find liblapack. Search for the correct paths"); + int jobvl; int jobvr; +/* + JOBVL: CHARACTER + JOBVR: CHARACTER + = 'N': left eigenvectors of A are not computed; + = 'V': left eigenvectors of A are computed. +*/ + int m,n,lda,ldvl,ldvr,lwork,info; - + var dgeev = dlfunc(lapack, "dgeev_", "vpppppppppppppp"); if(dgeev == null || !dgeev) die("Not DGEEV"); var WORK=new 'NumMatrix2D'; jobvl=jobvr= ord('V', 0); + m=A.rows; - n=A.cols; - - if(m==n) - { + n=A.cols; + lda=max(1,n); WR.resize(n,n); - WI.resize(n,n); + WI.resize(n,n); ldvl=ldvr=max(1,n); VL.resize(ldvl,n); VR.resize(ldvr,n); lwork=max(1,n*4); WORK.resize(lwork,lwork); - } - + dgeev(jobvl,jobvr,n,A,lda,WR,WI,VL,ldvl,VR,ldvr,WORK,lwork,info); + +/* + INFO is INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value. + > 0: if INFO = i, the QR algorithm failed to compute all the + eigenvalues, and no eigenvectors have been computed; + elements i+1:N of WR and WI contain eigenvalues which + have converged. + */ + return info; } + +//computes the max value + function max(var a,var b) { return a>b?a:b; } } -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'NumMatrix2D'; - var wr=new 'NumMatrix2D'; - var wi=new 'NumMatrix2D'; - var vl=new 'NumMatrix2D'; - var vr=new 'NumMatrix2D'; - var info; - a.initialize_from_args(3, 3, - 11.0, 2.0, 3.0, - 4.0, 5.0, 6.0, - 7.0, 1.0, 2.0); - - using dgeev_func_eve.dgeev_exec_eve; - info=dgeev_exec_eve(a,wr,wi,vl,vr); - - if(info==0) - { - say("right eigenvectors=\n",vr); - say("left eigenvectors=\n",vl); - } -} - -*/
50 src/LAPACK/inverse.winxed
 @@ -1,16 +1,33 @@ +//DGETRI computes the inverse of a matrix using the LU factorization computed by DGETRF. +//This method inverts U and then computes inv(A) by solving the system +//inv(A)*L = inv(U) for inv(A). + + \$load "Double.pbc"; namespace dgetri_func{ const int PRINT_DEBUG_STUFF = 0; +//a is DOUBLE PRECISION array, dimension (LDA,N) + function dgetri_exec(var a) { var pla = loadlib("linalg_group"); var lapack = loadlib('liblapack.so'); if (lapack == null || !lapack) die("Cannot find liblapack"); + say("Given Matrix:"); + say(a); + +/* +A is DOUBLE PRECISION array, dimension (LDA,N) + On entry, the factors L and U from the factorization + A = P*L*U as computed by DGETRF. + On exit, if INFO = 0, the inverse of the original matrix A. +*/ + int n,m,ipiv_size; m=a.rows; n=a.cols; @@ -20,8 +37,6 @@ function dgetri_exec(var a) for(i=0;i
42 src/LAPACK/linearequation.winxed
 @@ -1,7 +1,21 @@ +//DSPSV computes the solution to a real system of linear equations +// A * X = B, +//where A is an N-by-N symmetric matrix stored in packed format and Xand B are N-by-NRHS matrices. +//where U (or L) is a product of permutation and unit upper (lower) triangular matrices, D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then used to solve the system of equations A * X = B. + + namespace dspsv_func{ const int PRINT_DEBUG_STUFF = 0; +/* + A: DOUBLE PRECISION array, dimension (LDA,N) + B: DOUBLE PRECISION array, dimension (LDB,NRHS) + UPLO: CHARACTER + 'U': Upper triangle of A is stored + 'L': Lower triangle of A is stored. +*/ + function dspsv_exec(var A,var b,int uplo) { @@ -83,31 +97,3 @@ function debug(var matrix, string msg, var args [slurpy]) } -/* - -function main[main](var args) -{ - var pla = loadlib("linalg_group"); - var a = new 'NumMatrix2D'; - var b=new 'NumMatrix2D'; - var ap=new 'NumMatrix2D'; - int uplo; - uplo=85; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - b.initialize_from_args(3, 1, - 1.0, - 1.0, - 1.0); - using dspsv_func.dspsv_exec; - ap=dspsv_exec(a,b,uplo); - - say("AP="); - say(ap); - say("B="); - say(b); -} - -*/
61 src/LAPACK/multiplication.winxed
 @@ -1,7 +1,14 @@ +//MULTIPLICATION performs matrix multiplication using LAPACK library dgemm + namespace multiplication{ const int PRINT_DEBUG_STUFF = 0; +/* +A : DOUBLE PRECISION array of DIMENSION ( LDA, ka ) +B : DOUBLE PRECISION array of DIMENSION ( LDB, kb ) +*/ + function matmul(var A,var B) { @@ -34,38 +41,9 @@ function matmul(var A,var B) beta.resize(ldc,n); dgemm(transa,transb,m,n,k,alpha,A,lda,B,ldb,beta,C,ldc); return C; +//C is DOUBLE PRECISION array of DIMENSION ( LDC, n ) } -/* - int m1,n1,m2,n2; - m1=A.rows; - n1=A.cols; - m2=B.rows; - n2=B.cols; - if(n1==m2) - { - var pla = loadlib("linalg_group"); - var R=new 'NumMatrix2D'; - R.resize(m1,n2); - int i,j,k; - for(i=0;i
49 src/LAPACK/orthoqrfactorisation.winxed
 @@ -1,7 +1,18 @@ +//DORGQR generates an M-by-N real matrix Q with orthonormal columns,which is defined as the first N columns of a product of K lementary + reflectors of order M + + namespace dorgqr_func{ const int PRINT_DEBUG_STUFF = 0; +/* +A : DOUBLE PRECISION array, dimension (LDA,N) +K : INTEGER + The number of elementary reflectors whose product defines the + matrix Q. N >= K >= 0. +*/ + function dorgqr_exec(var A,int k) { @@ -47,18 +58,35 @@ function dorgqr_exec(var A,int k) for(i=0;ib?a:b; } +//computes the min value + function min(var a,var b) { return a>b?b:a; @@ -73,24 +101,3 @@ function debug(var matrix, string msg, var args [slurpy]) } } - -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'NumMatrix2D'; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - using dorgqr_func.dorgqr_exec; - int info; - info=dorgqr_exec(a,2); - if(info==0) - say(a); - say(""); -} - -*/
42 src/LAPACK/qrfactorisation.winxed
 @@ -1,7 +1,14 @@ +// DGEQRF computes a QR factorization of a real M-by-N matrix A: +// A = Q * R. + namespace dgeqrf_func{ const int PRINT_DEBUG_STUFF = 0; +/* +A : DOUBLE PRECISION array, dimension (LDA,N) +*/ + function dgeqrf_exec(var A) { @@ -32,11 +39,25 @@ function dgeqrf_exec(var A) var dgeqrf = dlfunc(lapack, "dgeqrf_", "vpppppppp"); if(dgeqrf == null || !dgeqrf) die("Not DGEQRF"); - + +/* + A is DOUBLE PRECISION array, dimension (LDA,N) + On entry, the M-by-N matrix A. + On exit, the elements on and above the diagonal of the array + contain the min(M,N)-by-N upper trapezoidal matrix R (R is + upper triangular if m >= n); the elements below the diagonal, + with the array TAU, represent the orthogonal matrix Q as a + product of min(m,n) elementary reflectors +*/ dgeqrf(m,n,A,lda,tau,work,lwork,info); return info; +/* + INFO is INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value +*/ } @@ -60,22 +81,3 @@ function debug(var matrix, string msg, var args [slurpy]) } -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'NumMatrix2D'; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - using dgeqrf_func.dgeqrf_exec; - int info; - info=dgeqrf_exec(a); - if(info==0) - say(a); -} - -*/
86 src/LAPACK/singlevaluedecomposition.winxed
 @@ -1,7 +1,33 @@ +/* +DGESVD computes the singular value decomposition (SVD) of a real M-by-N matrix A, optionally computing the left and/or right singular vectors. The SVD is written + + A = U * SIGMA * transpose(V) + + where SIGMA is an M-by-N matrix which is zero except for its min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and + V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A. +*/ + namespace dgesvd_func{ const int PRINT_DEBUG_STUFF = 0; +/* +A : DOUBLE PRECISION array, dimension (LDA,N) +JOBU : CHARACTER +JOBVT: CHARACTER + Specifies options for computing all or part of the matrix U: + = 'A': all M columns of U are returned in array U: + = 'S': the first min(m,n) columns of U (the left singular + vectors) are returned in the array U; + = 'O': the first min(m,n) columns of U (the left singular + vectors) are overwritten on the array A; + = 'N': no columns of U (no left singular vectors) are + computed. + S : DOUBLE PRECISION array, dimension (min(M,N)) + U : DOUBLE PRECISION array, dimension (LDU,UCOL) + VT : DOUBLE PRECISION array, dimension (LDVT,N) +*/ + function dgesvd_exec(var A,int jobu,int jobvt,var S,var U,var VT) { @@ -39,9 +65,33 @@ function dgesvd_exec(var A,int jobu,int jobvt,var S,var U,var VT) if(dgesvd == null || !dgesvd) die("Not DGESVD"); +/* + A is DOUBLE PRECISION array, dimension (LDA,N) + On entry, the M-by-N matrix A. + On exit, + if JOBU = 'O', A is overwritten with the first min(m,n) + columns of U (the left singular vectors, + stored columnwise); + if JOBVT = 'O', A is overwritten with the first min(m,n) + rows of V**T (the right singular vectors, + stored rowwise); + if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A + are destroyed. +*/ + dgesvd(jobu,jobvt,m,n,A,lda,S,U,ldu,VT,ldvt,work,lwork,info); return info; +/* + INFO is INTEGER + = 0: successful exit. + < 0: if INFO = -i, the i-th argument had an illegal value. + > 0: if DBDSQR did not converge, INFO specifies how many + superdiagonals of an intermediate bidiagonal form B + did not converge to zero. See the description of WORK + above for details. +*/ + } function max(var a,var b) @@ -63,39 +113,3 @@ function debug(var matrix, string msg, var args [slurpy]) } } - -/* - -function main[main](var args) -{ - - var pla = loadlib("linalg_group"); - var a = new 'NumMatrix2D'; - var S = new 'NumMatrix2D'; - var U = new 'NumMatrix2D'; - var VT = new 'NumMatrix2D'; - a.initialize_from_args(3, 3, - 1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, 0.0, 1.0); - - int jobu,jobvt; - jobu=jobvt= ord('A', 0); - - using dgesvd_func.dgesvd_exec; - int info; - - info=dgesvd_exec(a,jobu,jobvt,S,U,VT); - - if(info == 0) - { - say("successful"); - say(a); - say(S); - say(U); - say(VT); - } - -} - -*/

0 comments on commit `0baee2c`

Please sign in to comment.