/*-*-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
	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;.     */

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],SPECIAL))$

SSTATUS(FEATURE,EIGEN)$

HERMITIANMATRIX:FALSE$

NONDIAGONALIZABLE:FALSE$

KNOWNEIGVALS:FALSE$

KNOWNEIGVECTS:FALSE$

LISTEIGVECTS:[]$

LISTEIGVALS:[]$

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

COLUMNVECTOR(X):=TRANSPOSE(MATRIX(X))$


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


EIGENVALUES(MAT):=
		BLOCK([DIMNSN,LISTALL,SOLUTION,MULTIPLICITIES,SOLVEEXPLICIT,
	        dummy:?GENSYM(),INDEX2],
		LISTALL:[],
		SOLVEEXPLICIT:TRUE,
		DIMNSN:LENGTH(MAT),
		SOLUTION:BLOCK([PROGRAMMODE:TRUE],
	                       SOLVE(CHARPOLY(MAT,DUMMY),DUMMY)),
		IF SOLUTION=[] THEN 
		(PRINT("SOLVE is unable to find the roots of characteristic polynomial."),
		RETURN(LISTALL))
		ELSE (FOR INDEX2 THRU DIMNSN DO
		(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,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(NOTKNWN[INDEX1],UNKNOWNS),
		VECTR:COLUMNVECTOR(UNKNOWNS),MATRX:MAT.VECTR,NONDIAGONALIZABLE:FALSE,
		LISTALL:[EIGVALS],REALONLY:FALSE,
		FOR INDEX1 THRU COUNT DO 
		 (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 INDEX3 THRU %RNUM DO
		   (INTERN:SUBSTVECTK(%RNUM_LIST,INDEX3,INTERM),
                    UNIT:APPEND(UNIT,[INTERN])),

		  IF HERMITIANMATRIX AND %RNUM>1 THEN UNIT:GRAMSCHMIDT(UNIT),
		  LISTALL:APPEND(LISTALL,UNIT)),

		RETURN(LISTALL)))$


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):=
 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],
		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],
		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)$