--* From postmaster%watson.vnet.ibm.com@yktvmv.watson.ibm.com Thu Dec 1 18:00:46 1994 --* Received: from yktvmv-ob.watson.ibm.com by asharp.watson.ibm.com (AIX 3.2/UCB 5.64/930311) --* id AA20568; Thu, 1 Dec 1994 18:00:46 -0500 --* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3) --* with BSMTP id 9381; Thu, 01 Dec 94 18:00:52 EST --* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0" --* id --* for asbugs@watson; Thu, 01 Dec 94 18:00:51 -0500 --* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0" --* id 0981; Thu, 1 Dec 1994 18:00:48 EST --* Received: from pi.watson.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3) --* with TCP; Thu, 01 Dec 94 18:00:46 EST --* Received: by pi.watson.ibm.com (AIX 3.2/UCB 5.64/4.03) --* id AA15384; Thu, 1 Dec 1994 17:55:08 -0600 --* Date: Thu, 1 Dec 1994 17:55:08 -0600 --* From: rmc@pi.watson.ibm.com --* Message-Id: <9412012355.AA15384@pi.watson.ibm.com> --* To: asbugs@watson.ibm.com --* Subject: [2] Bug: 3 constraints not checked --@ Fixed by: SSD Mon Mar 20 08:24:33 EST 1995 --@ Tested by: none --@ Summary: Fixed type form creation for definitions with inferred types. -- Command line: asharp -Fx -lblas -llapack -lxlf -lasdem lazard.as adf.o -- Version: don't know -- Original bug file name: lazard.as -- ln -s $lapack/blas.a ./libblas.a -- ln -s $lapack/lapack.a ./liblapack.a -- lapack library on asharp is: /share/power/mathlib/lapack -- asharp -Fx -lblas -llapack -lxlf -lasdem lazard.as adf.o #include "aslib.as" #library demolib "asdem" #library svd "adf.aso" import from demolib; import from svd; DF ==> DoubleFloat; CDF ==> Complex DoubleFloat; SI ==> SingleInteger; MDF ==> Matrix_DoubleFloat; VDF ==> Vector_DoubleFloat; MCDF==> Matrix_Complex_DoubleFloat; VCDF==> Vector_Complex_DoubleFloat; RAD ==> PackedArray_DoubleFloat; RAC ==> PackedArray_PackedComplex_DoubleFloat; FtnInt ==> Record(ii: BSInt); FtnDouble ==> Record(aa: BDFlo); FtnChar ==> Record(bb: BChar); fortran (n: SI) : FtnInt == [n::BSInt]; fortran (x: DF) : FtnDouble == [x::BDFlo]; fortran (x: Character): FtnChar == [x::BChar]; asharp (n: FtnInt) : SingleInteger == n.ii::SingleInteger; import { -- LAPACK Singular Value Decomposition: A = u S vT. -- This updates A := v, S := S, B := uT B dgesvd: ( jobu : FtnChar, -- JOBU (input) CHARACTER*1 -- 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. -- jobvt : FtnChar, -- JOBVT (input) CHARACTER*1 -- Specifies options for computing all or part of the matrix -- V**T: -- = 'A': all N rows of V**T are returned in the array VT; -- = 'S': the first min(m,n) rows of V**T (the right singular -- vectors) are returned in the array VT; -- = 'O': the first min(m,n) rows of V**T (the right singular -- vectors) are overwritten on the array A; -- = 'N': no rows of V**T (no right singular vectors) are -- computed. -- -- JOBVT and JOBU cannot both be 'O'. -- m : FtnInt, -- M (input) INTEGER -- The number of rows of the input matrix A. M >= 0. -- n : FtnInt, -- N (input) INTEGER -- The number of columns of the input matrix A. N >= 0. -- a : RAD, -- A (input/output) 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. -- lda : FtnInt, -- LDA (input) INTEGER -- The leading dimension of the array A. LDA >= max(1,M). -- s : RAD, -- S (output) DOUBLE PRECISION array, dimension (min(M,N)) -- The singular values of A, sorted so that S(i) >= S(i+1). -- u : RAD, -- U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) -- (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. -- If JOBU = 'A', U contains the M-by-M orthogonal matrix U; -- if JOBU = 'S', U contains the first min(m,n) columns of U -- (the left singular vectors, stored columnwise); -- if JOBU = 'N' or 'O', U is not referenced. -- ldu : FtnInt, -- LDU (input) INTEGER -- The leading dimension of the array U. LDU >= 1; if -- JOBU = 'S' or 'A', LDU >= M. -- vt : RAD, -- VT (output) DOUBLE PRECISION array, dimension (LDVT,N) -- If JOBVT = 'A', VT contains the N-by-N orthogonal matrix -- V**T; -- if JOBVT = 'S', VT contains the first min(m,n) rows of -- V --*T (the right singular vectors, stored rowwise); -- if JOBVT = 'N' or 'O', VT is not referenced. -- ldvt : FtnInt, -- LDVT (input) INTEGER -- The leading dimension of the array VT. LDVT >= 1; if -- JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). -- work : RAD, -- WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) -- On exit, if INFO = 0, WORK(1) returns the optimal LWORK; -- if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged -- superdiagonal elements of an upper bidiagonal matrix B -- whose diagonal is in S (not necessarily sorted). B -- satisfies A = U * B * VT, so it has the same singular values -- as A, and singular vectors related by U and VT. -- lwork : FtnInt, -- LWORK (input) INTEGER -- The dimension of the array WORK. LWORK >= 1. -- LWORK >= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)-4). -- For good performance, LWORK should generally be larger. -- info : FtnInt -- INFO (output) 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. -- ) -> (); -- LAPACK Generalized eigenvalues: dgegv: ( jobvl : FtnChar, -- JOBVL (input) CHARACTER*1 -- = 'N': do not compute the left generalized eigenvectors; -- = 'V': compute the left generalized eigenvectors. -- jobvr : FtnChar, -- JOBVR (input) CHARACTER*1 -- = 'N': do not compute the right generalized eigenvectors; -- = 'V': compute the right generalized eigenvectors. -- n : FtnInt, -- N (input) INTEGER -- The number of rows and columns in the matrices A, B, VL, and -- VR. N >= 0. -- a : RAD, -- A (input/workspace) DOUBLE PRECISION array, dimension (LDA, N) -- On entry, the first of the pair of matrices whose -- generalized eigenvalues and (optionally) generalized -- eigenvectors are to be computed. -- On exit, the contents will have been destroyed. (For a -- description of the contents of A on exit, see "Further -- Details", below.) -- lda : FtnInt, -- LDA (input) INTEGER -- The leading dimension of A. LDA >= max(1,N). -- b : RAD, -- B (input/workspace) DOUBLE PRECISION array, dimension (LDB, N) -- On entry, the second of the pair of matrices whose -- generalized eigenvalues and (optionally) generalized -- eigenvectors are to be computed. -- On exit, the contents will have been destroyed. (For a -- description of the contents of B on exit, see "Further -- Details", below.) -- ldb : FtnInt, -- LDB (input) INTEGER -- The leading dimension of B. LDB >= max(1,N). -- alphar : RAD, -- ALPHAR (output) DOUBLE PRECISION array, dimension (N) alphai : RAD, -- ALPHAI (output) DOUBLE PRECISION array, dimension (N) beta : RAD, -- BETA (output) DOUBLE PRECISION array, dimension (N) -- -- On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will -- be the generalized eigenvalues. If ALPHAI(j) is zero, then -- the j-th eigenvalue is real; if positive, then the j-th and -- (j+1)-st eigenvalues are a complex conjugate pair, with -- ALPHAI(j+1) negative. -- -- Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) -- may easily over- or underflow, and BETA(j) may even be zero. -- Thus, the user should avoid naively computing the ratio -- alpha/beta. However, ALPHAR and ALPHAI will be always less -- than and usually comparable with norm(A) in magnitude, and -- BETA always less than and usually comparable with norm(B). -- vl : RAD, -- VL (output) DOUBLE PRECISION array, dimension (LDVL,N) -- If JOBVL = 'V', the left generalized eigenvectors. (See -- "Purpose", above.) Real eigenvectors take one column, -- complex take two columns, the first for the real part and -- the second for the imaginary part. Complex eigenvectors -- correspond to an eigenvalue with positive imaginary part. -- Each eigenvector will be scaled so the largest component -- will have abs(real part) + abs(imag. part) = 1, *except* -- that for eigenvalues with alpha=beta=0, a zero vector will -- be returned as the corresponding eigenvector. -- Not referenced if JOBVL = 'N'. -- ldvl : FtnInt, -- LDVL (input) INTEGER -- The leading dimension of the matrix VL. LDVL >= 1, and -- if JOBVL = 'V', LDVL >= N. -- vr : RAD, -- VR (output) DOUBLE PRECISION array, dimension (LDVR,N) -- If JOBVL = 'V', the right generalized eigenvectors. (See -- "Purpose", above.) Real eigenvectors take one column, -- complex take two columns, the first for the real part and -- the second for the imaginary part. Complex eigenvectors -- correspond to an eigenvalue with positive imaginary part. -- Each eigenvector will be scaled so the largest component -- will have abs(real part) + abs(imag. part) = 1, *except* -- that for eigenvalues with alpha=beta=0, a zero vector will -- be returned as the corresponding eigenvector. -- Not referenced if JOBVR = 'N'. -- ldvr : FtnInt, -- LDVR (input) INTEGER -- The leading dimension of the matrix VR. LDVR >= 1, and -- if JOBVR = 'V', LDVR >= N. -- work : RAD, -- WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) -- On exit, if INFO = 0, WORK(1) returns the optimal LWORK. -- lwork : FtnInt, -- LWORK (input) INTEGER -- The dimension of the array WORK. LWORK >= max(1,8*N). -- For good performance, LWORK must generally be larger. -- To compute the optimal value of LWORK, call ILAENV to get -- blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute: -- NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR; -- The optimal LWORK is: -- 2*N + MAX( 6*N, N*(NB+1) ). -- info : FtnInt -- INFO (output) INTEGER -- = 0: successful exit -- < 0: if INFO = -i, the i-th argument had an illegal value. -- = 1,...,N: -- The QZ iteration failed. No eigenvectors have been -- calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) -- should be correct for j=INFO+1,...,N. -- > N: errors that usually indicate LAPACK problems: -- =N+1: error return from DGGBAL -- =N+2: error return from DGEQRF -- =N+3: error return from DORMQR -- =N+4: error return from DORGQR -- =N+5: error return from DGGHRD -- =N+6: error return from DHGEQZ (other than failed -- iteration) -- =N+7: error return from DTGEVC -- =N+8: error return from DGGBAK (computing VL) -- =N+9: error return from DGGBAK (computing VR) -- =N+10: error return from DLASCL (various calls) ) -> () } from Foreign C; {Lazard(F:EuclideanDomain, nvars:SI):Exports == Implementation} where { Expon ==> HomogeneousDirectProduct nvars; DPoly ==> Polynomial(F, Expon); SI ==> SingleInteger; MF ==> Matrix F; Exports ==> with { lazardMat: (List DPoly) -> MF; ++ computes Lazard's matrix from a list of polys. lazardMat1: (DPoly, SI, List Expon) -> MF; ++ lazardMat1(poly, deg, monomList) computes section of Lazard's ++ matrix associated to poly in degree deg according to monomial ++ basis specified by monomList. computeDegree: (List DPoly) -> SI; ++ computes big degree for lazard (sum (di-1) + 1) monomials: (SI) -> List Expon; ++ monomials(deg) computes list of monomials in degrees up thru deg. } Implementation ==> add { import from DPoly; import from F; import from SingleInteger; monomials1(d:SI, n:SI):List List SI == { import from SI; d = 0 => list ([0 for i in 1..n]); n = 1 => list ([d]); ris:List List SI := empty(); for i in 0..d repeat for monom in monomials1(d - i,n - 1) repeat ris := cons(cons(i,monom),ris); ris; } directProduct(l:List SI):Expon == { v:Array SI := array(s for s in l); v pretend Expon -- needs a better constructor } monomials2(d:SI):List Expon == { import from List List SI; [directProduct(v) for v in monomials1(d,nvars)]; } monomials(deg:SI):List Expon == { import from List List Expon; reduce(concat,[monomials2(dd) for dd in 0..deg],[]) } position(m:Expon, l:List Expon):SI == { for i in 1.. for x in l repeat m=x => return i; return 0 } lazardMat1 (pol: DPoly, di:SI, ddmonomials:List Expon):MF == { import from List SI; import from List Expon; nr := #ddmonomials; lmonom := monomials(di); nc := #lmonom; cmat:MF := zero(nr,nc); while pol ~= 0 repeat { deg := degree pol; lc := leadingCoefficient pol; pol := reductum pol; for j in 1..nc for exp in lmonom repeat { (j1 := position(exp + deg, ddmonomials)) > 0 => cmat(j1,j) := lc } } cmat; } computeDegree(lpol :List DPoly):SI == { import from List SI; import from Expon; n1 := nvars + 1; ldeg := [sum degree f for f in lpol]; import from ListSort SI; ordeg := sort(>,cons(1,ldeg)); dd := (reduce(+,[deg for deg in ordeg for i in 1..n1],0)) - nvars; } lazardMat (lpol:List DPoly):MF == { #lpol < nvars => error "too many variables"; dd := computeDegree lpol; import from List Expon; ddmonomials := monomials(dd); import from List MF; reduce(horizConcat, [lazardMat1(pol, dd - sum degree pol, ddmonomials) for pol in lpol], zero(#ddmonomials,0)) } } } +++ following converts an asharp matrix to fortran style array +++ nrows and ncols give shape of the fortran array +++ row and col give placement positioning within the array (1 indexed) fortranMatrix(m:Matrix DF, nrows:SI, ncols:SI, fm:MDF(nrows, ncols), row:SI, col:SI):MDF(nrows,ncols) == { row := row-1; col := col-1; nr := nRows m; nc := nCols m; for i in 1..nr repeat for j in 1..nc repeat fm(row+i, col+j):= m(i,j); fm } fortranMatrix(m:Matrix DF,nrows:SI, ncols:SI):MDF(nrows, ncols) == { fortranMatrix(m, nrows, ncols, new(0), 1, 1) } subMatrix(nrows:SI, ncols:SI, mf:MDF(nrows, ncols), trow:SI, brow:SI, lcol:SI, rcol:SI, outrows:SI, outcols:SI): MDF(outrows, outcols) == { nmf :MDF(outrows, outcols) == new(0); for i in trow .. brow for ii in 1.. repeat for j in lcol..rcol for jj in 1.. repeat nmf(ii,jj) := mf(i,j); nmf } copy(nrows:SI, ncols:SI, mf:MDF(nrows, ncols)):MDF(nrows, ncols) == { tmf : MDF(nrows, ncols) := new(0); for i in 1..nrows repeat for j in 1..ncols repeat tmf(i,j) := mf(i,j); tmf } transpose(nrows:SI, ncols:SI, mf:MDF(nrows, ncols)):MDF(ncols, nrows) == { tmf : MDF(ncols, nrows) := new(0); for i in 1..nrows repeat for j in 1..ncols repeat tmf(j,i) := mf(i,j); tmf } asharpMatrix(fm:MDF(nrows, ncols), nrows:SI, ncols:SI, row:SI, col:SI, nrows2:SI, ncols2:SI):Matrix DF == { m:Matrix DF := zero(nrows2, ncols2); row := row-1; col := col-1; for j in 1..ncols2 repeat { for i in 1..nrows2 repeat { m(i,j) := fm(i+row, j+col) } } m } matrixProduct( arows:SI, acols:SI, bcols:SI, a:MDF(arows,acols), b:MDF(acols,bcols)):MDF(arows,bcols) == { ans : MDF(arows,bcols) := new 0; for i in 1..arows repeat { for j in 1..bcols repeat { for k in 1..acols repeat { ans(i,j) := ans(i,j) + a(i,k)*b(k,j); } } } ans } innerProducts(nrows:SI, lv:MCDF(nrows, nrows), mf:MDF(nrows, nrows), rv:MCDF(nrows,nrows)): VCDF(nrows) == { v:VCDF(nrows):=new(0); tempvec:VCDF(nrows) := new(0); for k in 1..nrows repeat { for i in 1..nrows repeat { tempvec(i) := 0; for j in 1..nrows repeat { tempvec(i) := tempvec(i) + (mf(i,j)::Complex DF)*rv(j,k); } } v(k) := 0; for ii in 1..nrows repeat v(k) := v(k) + tempvec(ii)*lv(ii,k); } v } R ==> DoubleFloat; hdp ==> HomogeneousDirectProduct; differentiate(nvars:SI, R:Field, p:Polynomial(R, hdp nvars), varind:SI): Polynomial(R, hdp nvars) == { zero? p => 0; import from hdp nvars; varpol := var unitVector varind; qr := monicDivide(p, varpol); poly ==> Polynomial(R, hdp nvars); import from Record(quotient:poly, remainder:poly); differentiate(nvars, R, qr.quotient, varind)*varpol + qr.quotient } import from SingleInteger; singularPts(p:Polynomial(DF, hdp 2), eps:DF):List List CDF == { lp : List Polynomial(DF, hdp 2) := [p, differentiate(2, DF, p, 1), differentiate(2, DF, p, 2)]; solve(2, lp, eps) } solve(nvars:SI, lp:List Polynomial(DF, hdp nvars), eps :DF): List List CDF == { poly ==> Polynomial(DF, hdp nvars); import from hdp nvars, poly; R ==> DF; -- create numerical matrix if fortran form "mf" -- import from Lazard(DF, nvars); m:Matrix R := lazardMat lp; nrows:SI == nRows m; ncols:SI == nCols m; print << "Polynomial system is: " << newline; for ll in lp repeat print << ll << newline; print << newline; import from List hdp nvars; -- ldm == max(nrows, ncols); ldm == nrows; mf : MDF(ldm, ncols) := new(0); mf := fortranMatrix(m, ldm, ncols, mf, 1, 1); -- create adjoined ui matrix b and call svd dd :SI := computeDegree lp; lmon := monomials dd; lvar :List poly := cons(1, [var unitVector i for i in 1..nvars]); nv == #lvar; blist : List Matrix DF := [lazardMat1(v, dd-1, lmon) for v in lvar]; nc:SI == nCols first blist; nb:SI == nv * nc; -- ldb == max(nrows, nb); ldb == nrows; b : MDF(ldb, nb) := new(0); for bb in blist for i in 0.. repeat b:=fortranMatrix(bb, ldb, nb, b, 1, 1 + i*nc); s:VDF(ncols):=callSVD(nrows, ncols, ldm, mf, nb, ldb, b); print << "singularValues are: " << newline; print << s << newline; seps := s.1 * eps; rank:SI := ncols; for i in 1..ncols repeat { if s(i) < seps then { rank := i-1; break; } } r:SI == nrows - rank; print << "System seems to have " << r << " geometric solutions." << newline; -- reduce to rxr matrices using SVD bigulist: List MDF(r,nc) := [subMatrix(ldb, nb, b, nrows-r+1, nrows, 1+i*nc, i*nc + nc, r, nc) for i in 0..nvars]; -- print << "bigulist = " << newline << bigulist << newline; ranbigu:MDF(nc,r) := transpose(r,nc, randomCombine(r, nc, bigulist)); -- print << "ranbigu = " << newline << ranbigu << newline; uadjoin : MDF(nc, nv*r) := new 0; for u in bigulist for k in 0.. repeat for i in 1..r repeat for j in 1..nc repeat uadjoin(j, k*r+i) := u(i,j); -- print << "uadjoin = " << newline << uadjoin << newline; callSVD(nc, r, nc, ranbigu, nv*r, nc, uadjoin); ulist: List MDF(r,r) := [subMatrix(nc, nv*r, uadjoin, 1, r, 1 + i*r, i*r + r, r, r) for i in 0..nvars]; -- print << "ulist = " << ulist << newline; -- compute eigenvectors for generic ui's and return eigenvalues u1 := randomCombine(r, r, ulist); u2 := randomCombine(r, r, ulist); print << "u1 is " << newline << u1 << newline; print << "u2 is " << newline << u2 << newline; local zl,zr :MCDF(r,r); z := callGEigen(r, u1, u2); (zl, zr) := z; import from List VCDF(r); wlist: List List CDF := normalize(r, [innerProducts(r, zl, uu, zr) for uu in ulist]); print << "Normalized solutions are :" << newline; for ww in wlist repeat print << ww << newline; wlist } normalize(r:SI, w:List VCDF(r)):List List CDF == { print << newline; print << "Eigenvalues are: " << newline; for ww in w repeat print << ww << newline; print << newline; [(mc:=maxCoord(r,w,i); if mc=0 then mc:=1; [uu(i)/mc for uu in w]) for i in 1..r] } maxCoord(r:SI, wlist:List VCDF(r), i:SI):CDF == { maxNorm:DF := 0; maxCoord:CDF := 0; for w in wlist repeat if norm(w.i)>maxNorm then { maxCoord := w.i; maxNorm := norm maxCoord } maxCoord } randomCombine(r:SI, c:SI, ulist: List MDF(r,c)): MDF(r,c) == { import from RandomNumberSource; den:DF := size :: DF; ranu : MDF(r,c) := new 0; for uu in ulist repeat { coef := random() ::DF / den; for i in 1..r repeat for j in 1..c repeat ranu(i,j) := ranu(i,j) + coef * uu(i,j); } ranu } callSVD(m:SI, n:SI, lda:SI, a:MDF(lda,n), nb:SI, ldb:SI, b:MDF(ldb,nb)): VDF(n) == { data xx ==> xx pretend RAD; -- print << "Inside CallSVD " << newline << "m is " << m << newline; -- print << "n is " << n << newline << "a is " << newline << a << newline; -- print << "b is " << newline; import from Character; jobu : FtnChar := fortran char "A"; -- form the whole U matrix jobvt: FtnChar := fortran char "N"; -- don't form the V matrix s : VDF(n) := new 0; -- singular values vector ldu : SI == m; -- ldb must be >= m. u : MDF(ldu,m) := new 0; -- U matrix ldvt : SI == 1; vt : MDF(ldvt,1) := new 0; -- dummy V matrix lwork : SI == max(3*min(m,n) + max(m,n), 5*min(m,n) - 4); work : VDF(lwork) := new 0; info : FtnInt := fortran 0; dgesvd( jobu, jobvt, fortran m, fortran n, data a, fortran lda, data s, data u, fortran ldu, data vt, fortran ldvt, data work, fortran lwork, info); if not zero? asharp info then error "Error return from dgesvd, info = ", info; print << "Inside CallSVD, u is " << newline << u << newline; -- One advantage of ESSL is that it did the multiplication of U'*b for us. newb : MDF(m,nb) := new 0; for i in 1..m repeat { for j in 1..nb repeat newb(i,j) := b(i,j); } uk:MDF(m,ldu) := transpose(ldu,m,u); ut:MDF(ldb,m) := new 0; for i in 1..m repeat { for j in 1..m repeat ut(i,j) := uk(i,j); } -- Side effect on return... b := matrixProduct(ldb, m, nb, ut, newb); print << "Inside CallSVD, b is " << newline << b << newline; s } callGEigen(r: SI, a:MDF(r,r), b:MDF(r,r)):(MCDF(r,r), MCDF(r,r)) == { data xx ==> xx pretend RAD; import from Character; jobvl : FtnChar := fortran char "V"; -- Compute left eigenvectors jobvr : FtnChar := fortran char "V"; -- Compute right eigenvectors alpha: VCDF(r) := new 0; -- Generalized eigenpairs ((alphar+i*alphai), beta) alphar:VDF(r) := new 0; alphai:VDF(r) := new 0; beta: VDF(r) := new 0; vl: MDF(r,r) := new 0; -- Real matrices for left and right eigenvectors. vr: MDF(r,r) := new 0; -- LAPACK uses real arithmetic. zl: MCDF(r,r) := new 0; -- Complex matrices expected outside. zr: MCDF(r,r) := new 0; -- zl = left eigenvectors, zr = right eigenvectors. r8: SI == 8*r; -- Easy estimate of work array size. work: VDF(r8) := new 0; info : FtnInt := fortran 0; print << "a is " << newline << a << newline; print << "b is " << newline << b << newline; dgegv( jobvl, jobvr, fortran r, data a, fortran r, data b, fortran r, data alphar, data alphai, data beta, data vl, fortran r, data vr, fortran r, data work, fortran r8, info); -- check the error return info code if not zero? asharp info then error "Error return from dgegv, info = ", info ; -- convert alphar + i*alphai to alpha, and convert vl, vr to zl, zr. for i in 1..r repeat { alpha(i) := complex( alphar(i), alphai(i) ) } j : SI := 0; while j < r repeat { j := j + 1; {zero? alphai(j) => for i in 1..r repeat { zl(i,j) := complex(vl(i,j), 0); zr(i,j) := complex(vr(i,j), 0); } -- otherwise for i in 1..r repeat { zl(i,j) := complex(vl(i,j), vl(i,j+1)); zl(i,j+1) := complex(vl(i,j), -vl(i,j+1)); zr(i,j) := complex(vr(i,j), vr(i,j+1)); zr(i,j+1) := complex(vr(i,j), -vr(i,j+1)); } j := j+1; } } print << "vl is " << newline << vl << newline; print << "vr is " << newline << vr << newline; print << "zl is " << newline << zl << newline; print << "zr is " << newline << zr << newline; -- need to add special handling for multiple eigenvalues (zl, zr) } test1():() =={ import from SingleInteger; R ==> DoubleFloat; import from R; -- create list of polynomials ---- import from hdp 2; poly ==> Polynomial(R, hdp 2); import from poly; x:= var unitVector 1; y:= var unitVector 2; p1 := x*x + x*y+ 2.0*x + y -1; q1 := x*x - y*y + 3.0*x + 2.0*y -1; lp:List poly := [p1, q1]; -- solve and print solutions wlist : List List CDF := solve(2, lp, .000001) } test1(); test2():() =={ import from SingleInteger; R ==> DoubleFloat; import from R; -- create polynomials ---- import from hdp 2; poly ==> Polynomial(R, hdp 2); import from poly; x:= var unitVector 1; y:= var unitVector 2; p:= 4.0*y^4+17.0*x^2*y^2+1.307*x*y^2-19.572938*y^2+4.0*x^4+5.228*x^3 -18.29175*x^2-5.228*x+15.29175*1; -- solve and print solutions singularPts(p, .00001) } test2();