/* modified for DOE MACSYMA with define_variable */

/*	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
	TO WASTE TIME (AND/OR PAPER...) THE FASTLOADABLE VERSION IS ON THE FILE
	EIGEN FASL DSK:SHARE;.  YOU CAN LOAD THE LATTER USING MACSYMA'S 
	LOADFILE COMMAND, I.E. LOADFILE(EIGEN,FASL,DSK,SHARE);. THE FUNCTIONS
	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),
        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)))$

/* COMMENTED OUT OF DOE MACSYMA
SSTATUS(FEATURE,EIGEN)$

HERMITIANMATRIX:FALSE$

NONDIAGONALIZABLE:FALSE$

KNOWNEIGVALS:FALSE$

KNOWNEIGVECTS:FALSE$

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)))$

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)],
		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
		))),
		COUNT:COUNT+1,LISTALL:ENDCONS(UNIT,LISTALL)),
		RETURN(LISTALL)))$


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)),
		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."),
		RETURN(LISTALL))
		ELSE (FOR INDEX2 THRU DIMNSN DO
		(DIMNSN:MI(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(),
		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 COUNT DO 
		(COUNT:MI(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 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].
   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)))$


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)
		ELSE (LISTALL:[PART(LISTUEVEC,1)],
		FOR INDEX1:2 THRU LENGTH(LISTUEVEC) DO 
		(UNIT:PART(LISTUEVEC,INDEX1),
		UNIT:RATSIMP(UNITVECTOR(UNIT)),
		LISTALL:ENDCONS(UNIT,LISTALL)),
		RETURN(LISTALL)))$


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))
		ELSE LEFTMATRIX:RIGHTMATRIX^^-1,
		RETURN(LISTUEVEC)))$


CONJ(X):=CONJUGATE(X)$

INPROD(X,Y):=INNERPRODUCT(X,Y)$

UVECT(X):=UNITVECTOR(X)$

COVECT(X):=COLUMNVECTOR(X)$

GSCHMIT(X):=GRAMSCHMIDT(X)$

EIVALS(MAT):=EIGENVALUES(MAT)$

EIVECTS(MAT):=EIGENVECTORS(MAT)$

UEIVECTS(MAT):=UNITEIGENVECTORS(MAT)$

SIMTRAN(MAT):=SIMILARITYTRANSFORM(MAT)$