Does anybody mind if I remove share/matrix/eigeni.mac ?
eigeni.mac hasn't been changed in any substantial way since the
beginning of time (May 8, 2000). There is also eigen.mac which has been
changed in various ways. It appears that originally eigen.mac was
supposed to be compiled and eigeni.mac was supposed to be interpreted;
there are no functional differences between the two that I can see.
There are no documentation references to eigeni.mac that I can find.
I've attached a diff of eigen.mac and eigeni.mac as they were in the
ur-commit 9e83265775.
best
Robert Dodier
PS.
--- share/matrix/eigen.mac 2013-10-06 20:13:08.409082814 -0700
+++ share/matrix/eigeni.mac 2013-10-06 20:13:08.409082814 -0700
@@ -1,5 +1,4 @@
-/* modified for DOE MACSYMA with define_variable */
-
+/*-*-macsyma-*-*/
/* THIS IS THE FILE EIGEN > DSK:SHARE;.
THIS IS THE SOURCE CODE FOR THE NEW EIGEN PACKAGE AND IT IS
MACSYMA BATCHABLE, I.E. BATCH(EIGEN,>,DSK,SHARE);. IF YOU DO NOT WANT
@@ -9,20 +8,15 @@
ARE DESCRIBED IN THE FILE EIGEN USAGE DSK:SHARE;, AND THE DEMO FILE IN
WHICH THE FUNCTIONS ARE DEMONSTRATED IS EIGEN DEMO DSK:SHARE;. */
-/* commented out of DOE MACSYMA
EVAL_WHEN(TRANSLATE_FILE,
MODEDECLARE([HERMITIANMATRIX,NONDIAGONALIZABLE,KNOWNEIGVALS,
KNOWNEIGVECTS],BOOLEAN,
- [INDEX1,INDEX2,INDEX3,INDEX4,DIMNSN,COUNT,%RNUM],FIXNUM),
+ [INDEX1,INDEX2,INDEX3,INDEX4,DIMNSN,COUNT,
+ %RNUM],FIXNUM),
DECLARE([HERMITIANMATRIX,NONDIAGONALIZABLE,KNOWNEIGVALS,
- KNOWNEIGVECTS,LISTEIGVECTS,LISTEIGVALS,%RNUM,LISTARITH,PROGRAMMODE,
- ALGEBRAIC,REALONLY,MULTIPLICITIES,SOLVEEXPLICIT],SPECIAL))$ */
-
-EVAL_WHEN([TRANSLATE,BATCH,DEMO,LOAD,LOADFILE],
-MI(VAR)::=BUILDQ([VAR],MODE_IDENTITY(FIXNUM,VAR)),
-DV(VAR)::=BUILDQ([VAR],DEFINE_VARIABLE(VAR,FALSE,BOOLEAN)))$
+ KNOWNEIGVECTS,LISTEIGVECTS,LISTEIGVALS,%RNUM,LISTARITH,
+ PROGRAMMODE],SPECIAL))$
-/* COMMENTED OUT OF DOE MACSYMA
SSTATUS(FEATURE,EIGEN)$
HERMITIANMATRIX:FALSE$
@@ -36,40 +30,24 @@
LISTEIGVECTS:[]$
LISTEIGVALS:[]$
-*/
-
-DV(HERMITIANMATRIX)$ DV(NONDIAGONALIZABLE)$ DV(KNOWNEIGVALS)$
-DV(KNOWNEIGVECTS)$
-
-DEFINE_VARIABLE(LISTEIGVECTS,[],LIST)$
-DEFINE_VARIABLE(LISTEIGVALS,[],LIST)$
-DEFINE_VARIABLE(RIGHTMATRIX,[],ANY)$
-DEFINE_VARIABLE(LEFTMATRIX,[],ANY)$
CONJUGATE(X):=SUBLIS('[%I=-%I],X)$
INNERPRODUCT(X,Y):=BLOCK([LISTARITH],LISTARITH:TRUE,RATSIMP(CONJUGATE(X).Y))$
-/*
-UNITVECTOR(X):=BLOCK([LISTARITH,INTRN],LISTARITH:TRUE,INTRN:INNERPRODUCT(X,X),
-INTRN:SQRT(INTRN),X/INTRN)$
-*/
-UNITVECTOR(X):=BLOCK([LISTARITH],LISTARITH:TRUE,X/SQRT(INNERPRODUCT(X,X)))$
+
+UNITVECTOR(X):=BLOCK([LISTARITH,INTRN],LISTARITH:TRUE,INTRN:INNERPRODUCT(X,X),INTRN:SQRT(INTRN),X/INTRN)$
COLUMNVECTOR(X):=TRANSPOSE(MATRIX(X))$
GRAMSCHMIDT(X):=
- BLOCK([LISTARITH,DIMNSN,LISTALL,INTERN,COUNT,DENOM,UNIT,INDEX1,
- INDEX2],
- MODE_DECLARE([DIMNSN,COUNT,INDEX1,INDEX2],FIXNUM,
- [LISTALL],LIST,[INTERN,DENOM,UNIT],ANY),
- LISTARITH:TRUE,DIMNSN:MI(LENGTH(X)),LISTALL:[PART(X,1)],
+ BLOCK([LISTARITH,DIMNSN,LISTALL,INTERN,COUNT,DENOM,UNIT,INDEX1,INDEX2],
+ LISTARITH:TRUE,DIMNSN:LENGTH(X),LISTALL:[PART(X,1)],
COUNT:1,IF DIMNSN=1 THEN RETURN(X)
ELSE (FOR INDEX1:2 THRU DIMNSN DO
(UNIT:PART(X,INDEX1),FOR INDEX2 THRU COUNT DO
(INTERN:PART(LISTALL,INDEX2),DENOM:INNERPRODUCT(INTERN,INTERN),
- UNIT:FACTOR(RATSIMP(UNIT-INNERPRODUCT(INTERN,UNIT)*INTERN/DENOM
- ))),
+ UNIT:FACTOR(RATSIMP(UNIT-INNERPRODUCT(INTERN,UNIT)*INTERN/DENOM))),
COUNT:COUNT+1,LISTALL:ENDCONS(UNIT,LISTALL)),
RETURN(LISTALL)))$
@@ -77,76 +55,90 @@
EIGENVALUES(MAT):=
BLOCK([DIMNSN,LISTALL,SOLUTION,MULTIPLICITIES,SOLVEEXPLICIT,
DUMMY:?GENSYM(),INDEX2],
- MODE_DECLARE([DIMNSN,INDEX2],FIXNUM,[LISTALL,SOLUTION],LIST,
- [DUMMY],ANY),
LISTALL:[],
SOLVEEXPLICIT:TRUE,
- DIMNSN:MI(LENGTH(MAT)),
+ DIMNSN:LENGTH(MAT),
SOLUTION:BLOCK([PROGRAMMODE:TRUE],
SOLVE(CHARPOLY(MAT,DUMMY),DUMMY)),
IF SOLUTION=[] THEN
- (PRINT(" "),PRINT("SOLVE is unable to find the roots of"),
- PRINT("the characteristic polynomial."),
+ (PRINT("SOLVE is unable to find the roots of characteristic polynomial."),
RETURN(LISTALL))
ELSE (FOR INDEX2 THRU DIMNSN DO
- (DIMNSN:MI(DIMNSN-PART(MULTIPLICITIES,INDEX2)+1),
+ (DIMNSN:DIMNSN-PART(MULTIPLICITIES,INDEX2)+1,
LISTALL:ENDCONS(RHS(PART(SOLUTION,INDEX2)),LISTALL)),
LISTALL:ENDCONS(MULTIPLICITIES,[LISTALL]),
RETURN(LISTALL)))$
+
EIGENVECTORS(MAT):=
- BLOCK([EQUATIONS,UNKNOWNS,SOLUTION,LISTALL,EIGVALS,DIMNSN,
- COUNT,VECTR,INDEX3,INDEX4,INDEX2,INDEX1,MATRX,MMATRX,notknwn,
- UNIT,MULTIPLICITIES,%RNUM,REALONLY,ALGEBRAIC,INTERM,INTERN],
- MODE_DECLARE([EQUATIONS,UNKNOWNS,SOLUTION,LISTALL,EIGVALS,
- UNIT,INTERM],LIST,
- [DIMNSN,COUNT,INDEX3,INDEX4,INDEX2,INDEX1],FIXNUM,
- [VECTR,MATRX,MMATRX,INTERN,NOTKNWN],ANY),
- UNKNOWNS:[],DIMNSN:MI(LENGTH(MAT)),
- COUNT:MI(DIMNSN),notknwn:?gensym(),
+ BLOCK([EQUATIONS,UNKNOWNS,SOLUTION,LISTALL,EIGVALS,DIMNSN,COUNT,VECTR,
+ INDEX3,INDEX4,INDEX2,INDEX1,NOTKNWN,MATRX,MMATRX,UNIT,MULTIPLICITIES,
+ %RNUM,REALONLY,INTERM,INTERN],
+ LOCAL(NOTKNWN),UNKNOWNS:[],DIMNSN:LENGTH(MAT),COUNT:DIMNSN,
IF KNOWNEIGVALS THEN EIGVALS:LISTEIGVALS
ELSE EIGVALS:EIGENVALUES(MAT),
IF EIGVALS=[] THEN (NONDIAGONALIZABLE:TRUE,RETURN(EIGVALS))
ELSE (MULTIPLICITIES:PART(EIGVALS,2),
- FOR INDEX1 THRU DIMNSN DO
- UNKNOWNS:ENDCONS(concat(notknwn,index1),UNKNOWNS),
- VECTR:COLUMNVECTOR(UNKNOWNS),MATRX:MAT.VECTR,
- NONDIAGONALIZABLE:FALSE,
- LISTALL:[EIGVALS],REALONLY:FALSE,ALGEBRAIC:TRUE,
+ FOR INDEX1 THRU DIMNSN DO UNKNOWNS:ENDCONS(NOTKNWN[INDEX1],UNKNOWNS),
+ VECTR:COLUMNVECTOR(UNKNOWNS),MATRX:MAT.VECTR,NONDIAGONALIZABLE:FALSE,
+ LISTALL:[EIGVALS],REALONLY:FALSE,
FOR INDEX1 THRU COUNT DO
- (COUNT:MI(COUNT-PART(MULTIPLICITIES,INDEX1)+1),
+ (COUNT:COUNT-PART(MULTIPLICITIES,INDEX1)+1,
MMATRX:MATRX-PART(EIGVALS,1,INDEX1)*VECTR,
EQUATIONS:[],
- FOR INDEX2 THRU DIMNSN DO
- EQUATIONS:CONS(MMATRX[INDEX2,1],EQUATIONS),%RNUM:0,
- SOLUTION:ALGSYS(EQUATIONS,UNKNOWNS),
- INTERM:MAP('RHS,SOLUTION[1]),
- UNIT:[],IF %RNUM#PART(MULTIPLICITIES,INDEX1)
- THEN NONDIAGONALIZABLE:TRUE,
+ FOR INDEX2 THRU DIMNSN DO EQUATIONS:CONS(MMATRX[INDEX2,1],EQUATIONS),%RNUM:0,
+ SOLUTION:ALGSYS(EQUATIONS,UNKNOWNS),INTERM:MAP(RHS,SOLUTION[1]),
+ UNIT:[],IF %RNUM#PART(MULTIPLICITIES,INDEX1) THEN NONDIAGONALIZABLE:TRUE,
FOR INDEX3 THRU %RNUM DO
(INTERN:SUBSTVECTK(%RNUM_LIST,INDEX3,INTERM),
UNIT:APPEND(UNIT,[INTERN])),
- IF UNIT=[] THEN
- (PRINT(" "),PRINT("ALGSYS failure: The eigenvector(s) for the",
- INDEX1,"th eigenvalue will be missing.")),
+
IF HERMITIANMATRIX AND %RNUM>1 THEN UNIT:GRAMSCHMIDT(UNIT),
LISTALL:APPEND(LISTALL,UNIT)),
+
RETURN(LISTALL)))$
-/* The first arg is of the form [r1,r2,r3].
+mEIGENVECTORS(MAT):=
+ BLOCK([EIGVALS:if knowneigvals then listeigvals
+ else eigenvalues(mat),
+ DIMNSN:length(mat),
+ UNIT],
+ IF EIGVALS=[] THEN (NONDIAGONALIZABLE:TRUE,RETURN(EIGVALS)),
+ NONDIAGONALIZABLE:FALSE,
+ cons(eigvals,
+ apply('append,
+ apply('map,
+ cons(lambda([val,multiplicity],
+ unit:args(transpose(%meig(echelon(mat-val*ident(dimnsn)),
+ dimnsn))),
+ if multiplicity#length(unit) then nondiagonalizable:true,
+ if hermitianmatrix then gramschmidt(unit) else unit),
+ eigvals)))))$
+
+%meig(mat,dimnsn):=
+ if mat=matrix() then matrix()
+ else
+ if mat[1,1]=0 then addcol(ematrix(dimnsn,1,1,1,1),
+ (%meig(submatrix(dimnsn,mat,1),
+ dimnsn-1),
+ addrow(matrix(0*%%[1]),%%)))
+ else (%meig(submatrix(1,mat,1),dimnsn-1),
+ addrow(-rest(mat[1]).%%,%%))$
+
+/* The first arg is of the form [r1,r2,r3]
We want to construct [r1=0,r2=1,r3=0] for example. */
-SUBSTVECTK(L,N,EXP):=(mode_declare(l,list,n,fixnum,exp,any),
- BLOCK([SUB_LIST:[],J:0],mode_declare(sub_list,list,j,fixnum),
- FOR VAR IN L DO (mode_declare(var,any),
- J:J+1,SUB_LIST:CONS(VAR = IF J=N THEN 1 ELSE 0,SUB_LIST)),
- SUBLIS(SUB_LIST,EXP)))$
+SUBSTVECTK(L,N,EXP):=
+ BLOCK([SUB_LIST:[],J:0],
+ FOR VAR IN L DO
+ (J:J+1,
+ SUB_LIST:CONS(VAR = IF J=N THEN 1 ELSE 0,SUB_LIST)),
+ SUBLIS(SUB_LIST,EXP))$
UNITEIGENVECTORS(MAT):=
BLOCK([LISTUEVEC,LISTALL,INDEX1,UNIT],
- mode_declare([listuevec,listall],list,index1,fixnum,unit,any),
IF KNOWNEIGVECTS THEN LISTUEVEC:LISTEIGVECTS
ELSE LISTUEVEC:EIGENVECTORS(MAT),
IF LISTUEVEC=[] THEN RETURN(LISTUEVEC)
@@ -160,13 +152,11 @@
SIMILARITYTRANSFORM(MAT):=
BLOCK([LISTVEC,LISTUEVEC],
- mode_declare([listvec,listuevec],list),
LISTUEVEC:UNITEIGENVECTORS(MAT),
IF NONDIAGONALIZABLE THEN RETURN(LISTUEVEC)
ELSE (LISTVEC:DELETE(PART(LISTUEVEC,1),LISTUEVEC),
- RIGHTMATRIX:TRANSPOSE(APPLY('MATRIX,LISTVEC)),
- IF HERMITIANMATRIX THEN
- LEFTMATRIX:CONJUGATE(TRANSPOSE(RIGHTMATRIX))
+ RIGHTMATRIX:TRANSPOSE(APPLY(MATRIX,LISTVEC)),
+ IF HERMITIANMATRIX THEN LEFTMATRIX:CONJUGATE(TRANSPOSE(RIGHTMATRIX))
ELSE LEFTMATRIX:RIGHTMATRIX^^-1,
RETURN(LISTUEVEC)))$