Author: neteler
Date: 2008-01-24 03:50:35 -0500 (Thu, 24 Jan 2008)
New Revision: 29816
Added:
grass-addons/m.eigensystem/
grass-addons/m.eigensystem/Makefile
grass-addons/m.eigensystem/README
grass-addons/m.eigensystem/awk.normal
grass-addons/m.eigensystem/awk.scale
grass-addons/m.eigensystem/balanc.f
grass-addons/m.eigensystem/balbak.f
grass-addons/m.eigensystem/cdiv.f
grass-addons/m.eigensystem/corr.7x7
grass-addons/m.eigensystem/elmhes.f
grass-addons/m.eigensystem/eltran.f
grass-addons/m.eigensystem/hqr.f
grass-addons/m.eigensystem/hqr2.f
grass-addons/m.eigensystem/m.eigensystem.html
grass-addons/m.eigensystem/main.f
grass-addons/m.eigensystem/matrix.7x7
grass-addons/m.eigensystem/rg.f
Log:
pseudo-ported m.eigensystem from GRASS 5
Added: grass-addons/m.eigensystem/Makefile
===================================================================
--- grass-addons/m.eigensystem/Makefile (rev 0)
+++ grass-addons/m.eigensystem/Makefile 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,25 @@
+MODULE_TOPDIR = ../..
+
+PGM = m.eigensystem
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: $(BIN)/$(PGM) html
+
+OBJ = main.o\
+ balanc.o \
+ balbak.o \
+ cdiv.o \
+ elmhes.o \
+ eltran.o \
+ hqr.o \
+ hqr2.o \
+ rg.o
+
+$(BIN)/$(PGM): $(OBJ)
+ $(FC) -o $@ $(OBJ)
+
+html:
+ $(MKDIR) $(GISBASE)/docs/html
+ $(INSTALL_DATA) $(PGM).html $(GISBASE)/docs/html/
+
Property changes on: grass-addons/m.eigensystem/Makefile
___________________________________________________________________
Name: svn:eol-style
+ native
Added: grass-addons/m.eigensystem/README
===================================================================
--- grass-addons/m.eigensystem/README (rev 0)
+++ grass-addons/m.eigensystem/README 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1 @@
+This code requires a fortran compiler
Added: grass-addons/m.eigensystem/awk.normal
===================================================================
--- grass-addons/m.eigensystem/awk.normal (rev 0)
+++ grass-addons/m.eigensystem/awk.normal 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,3 @@
+/E/ {i=0;sum=0.0;print}
+/V/ {i++;V[i]=$2;sum += $2*$2;print;n=0}
+/N/ {n++; print $0, V[n] / sqrt(sum), sqrt(sum)}
Added: grass-addons/m.eigensystem/awk.scale
===================================================================
--- grass-addons/m.eigensystem/awk.scale (rev 0)
+++ grass-addons/m.eigensystem/awk.scale 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,4 @@
+/V/ {print $0 " " $2/E; next}
+/N/ {print $0 " " $2/E; next}
+/E/ {E=sqrt($2)}
+{print}
Added: grass-addons/m.eigensystem/balanc.f
===================================================================
--- grass-addons/m.eigensystem/balanc.f (rev 0)
+++ grass-addons/m.eigensystem/balanc.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,166 @@
+ SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE) EIS0
+C EIS0
+ INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC EIS0
+ DOUBLE PRECISION A(NM,N),SCALE(N) EIS0
+ DOUBLE PRECISION C,F,G,R,S,B2,RADIX EIS0
+ LOGICAL NOCONV EIS0
+C EIS0
+C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, EIS0
+C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. EIS0
+C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). EIS0
+C EIS0
+C THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES EIS0
+C EIGENVALUES WHENEVER POSSIBLE. EIS0
+C EIS0
+C ON INPUT EIS0
+C EIS0
+C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL EIS0
+C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM EIS0
+C DIMENSION STATEMENT. EIS0
+C EIS0
+C N IS THE ORDER OF THE MATRIX. EIS0
+C EIS0
+C A CONTAINS THE INPUT MATRIX TO BE BALANCED. EIS0
+C EIS0
+C ON OUTPUT EIS0
+C EIS0
+C A CONTAINS THE BALANCED MATRIX. EIS0
+C EIS0
+C LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J) EIS0
+C IS EQUAL TO ZERO IF EIS0
+C (1) I IS GREATER THAN J AND EIS0
+C (2) J=1,...,LOW-1 OR I=IGH+1,...,N. EIS0
+C EIS0
+C SCALE CONTAINS INFORMATION DETERMINING THE EIS0
+C PERMUTATIONS AND SCALING FACTORS USED. EIS0
+C EIS0
+C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH EIS0
+C HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED EIS0
+C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS EIS0
+C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN EIS0
+C SCALE(J) = P(J), FOR J = 1,...,LOW-1 EIS0
+C = D(J,J), J = LOW,...,IGH EIS0
+C = P(J) J = IGH+1,...,N. EIS0
+C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, EIS0
+C THEN 1 TO LOW-1. EIS0
+C EIS0
+C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. EIS0
+C EIS0
+C THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN EIS0
+C BALANC IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS EIS0
+C K,L HAVE BEEN REVERSED.) EIS0
+C EIS0
+C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, EIS0
+C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY EIS0
+C EIS0
+C THIS VERSION DATED AUGUST 1983. EIS0
+C EIS0
+C ------------------------------------------------------------------EIS0
+C EIS0
+ RADIX = 16.0D0 EIS0
+C EIS0
+ B2 = RADIX * RADIX EIS0
+ K = 1 EIS0
+ L = N EIS0
+ GO TO 100 EIS0
+C .......... IN-LINE PROCEDURE FOR ROW AND EIS0
+C COLUMN EXCHANGE .......... EIS0
+ 20 SCALE(M) = J EIS0
+ IF (J .EQ. M) GO TO 50 EIS0
+C EIS0
+ DO 30 I = 1, L EIS0
+ F = A(I,J) EIS0
+ A(I,J) = A(I,M) EIS0
+ A(I,M) = F EIS0
+ 30 CONTINUE EIS0
+C EIS0
+ DO 40 I = K, N EIS0
+ F = A(J,I) EIS0
+ A(J,I) = A(M,I) EIS0
+ A(M,I) = F EIS0
+ 40 CONTINUE EIS0
+C EIS0
+ 50 GO TO (80,130), IEXC EIS0
+C .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE EIS0
+C AND PUSH THEM DOWN .......... EIS0
+ 80 IF (L .EQ. 1) GO TO 280 EIS0
+ L = L - 1 EIS0
+C .......... FOR J=L STEP -1 UNTIL 1 DO -- .......... EIS0
+ 100 DO 120 JJ = 1, L EIS0
+ J = L + 1 - JJ EIS0
+C EIS0
+ DO 110 I = 1, L EIS0
+ IF (I .EQ. J) GO TO 110 EIS0
+ IF (A(J,I) .NE. 0.0D0) GO TO 120 EIS0
+ 110 CONTINUE EIS0
+C EIS0
+ M = L EIS0
+ IEXC = 1 EIS0
+ GO TO 20 EIS0
+ 120 CONTINUE EIS0
+C EIS0
+ GO TO 140 EIS0
+C .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE EIS0
+C AND PUSH THEM LEFT .......... EIS0
+ 130 K = K + 1 EIS0
+C EIS0
+ 140 DO 170 J = K, L EIS0
+C EIS0
+ DO 150 I = K, L EIS0
+ IF (I .EQ. J) GO TO 150 EIS0
+ IF (A(I,J) .NE. 0.0D0) GO TO 170 EIS0
+ 150 CONTINUE EIS0
+C EIS0
+ M = K EIS0
+ IEXC = 2 EIS0
+ GO TO 20 EIS0
+ 170 CONTINUE EIS0
+C .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L .......... EIS0
+ DO 180 I = K, L EIS0
+ 180 SCALE(I) = 1.0D0 EIS0
+C .......... ITERATIVE LOOP FOR NORM REDUCTION .......... EIS0
+ 190 NOCONV = .FALSE. EIS0
+C EIS0
+ DO 270 I = K, L EIS0
+ C = 0.0D0 EIS0
+ R = 0.0D0 EIS0
+C EIS0
+ DO 200 J = K, L EIS0
+ IF (J .EQ. I) GO TO 200 EIS0
+ C = C + DABS(A(J,I)) EIS0
+ R = R + DABS(A(I,J)) EIS0
+ 200 CONTINUE EIS0
+C .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW .......... EIS0
+ IF (C .EQ. 0.0D0 .OR. R .EQ. 0.0D0) GO TO 270 EIS0
+ G = R / RADIX EIS0
+ F = 1.0D0 EIS0
+ S = C + R EIS0
+ 210 IF (C .GE. G) GO TO 220 EIS0
+ F = F * RADIX EIS0
+ C = C * B2 EIS0
+ GO TO 210 EIS0
+ 220 G = R * RADIX EIS0
+ 230 IF (C .LT. G) GO TO 240 EIS0
+ F = F / RADIX EIS0
+ C = C / B2 EIS0
+ GO TO 230 EIS0
+C .......... NOW BALANCE .......... EIS0
+ 240 IF ((C + R) / F .GE. 0.95D0 * S) GO TO 270 EIS0
+ G = 1.0D0 / F EIS0
+ SCALE(I) = SCALE(I) * F EIS0
+ NOCONV = .TRUE. EIS0
+C EIS0
+ DO 250 J = K, N EIS0
+ 250 A(I,J) = A(I,J) * G EIS0
+C EIS0
+ DO 260 J = 1, L EIS0
+ 260 A(J,I) = A(J,I) * F EIS0
+C EIS0
+ 270 CONTINUE EIS0
+C EIS0
+ IF (NOCONV) GO TO 190 EIS0
+C EIS0
+ 280 LOW = K EIS0
+ IGH = L EIS0
+ RETURN EIS0
+ END EIS0
Added: grass-addons/m.eigensystem/balbak.f
===================================================================
--- grass-addons/m.eigensystem/balbak.f (rev 0)
+++ grass-addons/m.eigensystem/balbak.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,75 @@
+ SUBROUTINE BALBAK(NM,N,LOW,IGH,SCALE,M,Z) EIS0
+C EIS0
+ INTEGER I,J,K,M,N,II,NM,IGH,LOW EIS0
+ DOUBLE PRECISION SCALE(N),Z(NM,M) EIS0
+ DOUBLE PRECISION S EIS0
+C EIS0
+C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK, EIS0
+C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. EIS0
+C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). EIS0
+C EIS0
+C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL EIS0
+C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING EIS0
+C BALANCED MATRIX DETERMINED BY BALANC. EIS0
+C EIS0
+C ON INPUT EIS0
+C EIS0
+C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL EIS0
+C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM EIS0
+C DIMENSION STATEMENT. EIS0
+C EIS0
+C N IS THE ORDER OF THE MATRIX. EIS0
+C EIS0
+C LOW AND IGH ARE INTEGERS DETERMINED BY BALANC. EIS0
+C EIS0
+C SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS EIS0
+C AND SCALING FACTORS USED BY BALANC. EIS0
+C EIS0
+C M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. EIS0
+C EIS0
+C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- EIS0
+C VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. EIS0
+C EIS0
+C ON OUTPUT EIS0
+C EIS0
+C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIS0
+C TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. EIS0
+C EIS0
+C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, EIS0
+C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY EIS0
+C EIS0
+C THIS VERSION DATED AUGUST 1983. EIS0
+C EIS0
+C ------------------------------------------------------------------EIS0
+C EIS0
+ IF (M .EQ. 0) GO TO 200 EIS0
+ IF (IGH .EQ. LOW) GO TO 120 EIS0
+C EIS0
+ DO 110 I = LOW, IGH EIS0
+ S = SCALE(I) EIS0
+C .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED EIS0
+C IF THE FOREGOING STATEMENT IS REPLACED BY EIS0
+C S=1.0D0/SCALE(I). .......... EIS0
+ DO 100 J = 1, M EIS0
+ 100 Z(I,J) = Z(I,J) * S EIS0
+C EIS0
+ 110 CONTINUE EIS0
+C ......... FOR I=LOW-1 STEP -1 UNTIL 1, EIS0
+C IGH+1 STEP 1 UNTIL N DO -- .......... EIS0
+ 120 DO 140 II = 1, N EIS0
+ I = II EIS0
+ IF (I .GE. LOW .AND. I .LE. IGH) GO TO 140 EIS0
+ IF (I .LT. LOW) I = LOW - II EIS0
+ K = SCALE(I) EIS0
+ IF (K .EQ. I) GO TO 140 EIS0
+C EIS0
+ DO 130 J = 1, M EIS0
+ S = Z(I,J) EIS0
+ Z(I,J) = Z(K,J) EIS0
+ Z(K,J) = S EIS0
+ 130 CONTINUE EIS0
+C EIS0
+ 140 CONTINUE EIS0
+C EIS0
+ 200 RETURN EIS0
+ END EIS0
Added: grass-addons/m.eigensystem/cdiv.f
===================================================================
--- grass-addons/m.eigensystem/cdiv.f (rev 0)
+++ grass-addons/m.eigensystem/cdiv.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,16 @@
+ SUBROUTINE CDIV(AR,AI,BR,BI,CR,CI) EIS0
+ DOUBLE PRECISION AR,AI,BR,BI,CR,CI EIS0
+C EIS0
+C COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI) EIS0
+C EIS0
+ DOUBLE PRECISION S,ARS,AIS,BRS,BIS EIS0
+ S = DABS(BR) + DABS(BI) EIS0
+ ARS = AR/S EIS0
+ AIS = AI/S EIS0
+ BRS = BR/S EIS0
+ BIS = BI/S EIS0
+ S = BRS**2 + BIS**2 EIS0
+ CR = (ARS*BRS + AIS*BIS)/S EIS0
+ CI = (AIS*BRS - ARS*BIS)/S EIS0
+ RETURN EIS0
+ END EIS0
Added: grass-addons/m.eigensystem/corr.7x7
===================================================================
--- grass-addons/m.eigensystem/corr.7x7 (rev 0)
+++ grass-addons/m.eigensystem/corr.7x7 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,8 @@
+7
+1.000 0.954 0.932 0.729 0.843 0.365 0.840
+0.954 1.000 0.952 0.795 0.877 0.390 0.862
+0.932 0.952 1.000 0.697 0.892 0.416 0.908
+0.729 0.795 0.697 1.000 0.819 0.459 0.673
+0.843 0.877 0.892 0.819 1.000 0.582 0.948
+0.365 0.390 0.416 0.459 0.582 1.000 0.552
+0.840 0.862 0.908 0.673 0.948 0.552 1.000
Added: grass-addons/m.eigensystem/elmhes.f
===================================================================
--- grass-addons/m.eigensystem/elmhes.f (rev 0)
+++ grass-addons/m.eigensystem/elmhes.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,98 @@
+ SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT) EIS3
+C EIS3
+ INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1 EIS3
+ DOUBLE PRECISION A(NM,N) EIS3
+ DOUBLE PRECISION X,Y EIS3
+ INTEGER INT(IGH) EIS3
+C EIS3
+C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES, EIS3
+C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. EIS3
+C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). EIS3
+C EIS3
+C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE EIS3
+C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS EIS3
+C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY EIS3
+C STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. EIS3
+C EIS3
+C ON INPUT EIS3
+C EIS3
+C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL EIS3
+C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM EIS3
+C DIMENSION STATEMENT. EIS3
+C EIS3
+C N IS THE ORDER OF THE MATRIX. EIS3
+C EIS3
+C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING EIS3
+C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, EIS3
+C SET LOW=1, IGH=N. EIS3
+C EIS3
+C A CONTAINS THE INPUT MATRIX. EIS3
+C EIS3
+C ON OUTPUT EIS3
+C EIS3
+C A CONTAINS THE HESSENBERG MATRIX. THE MULTIPLIERS EIS3
+C WHICH WERE USED IN THE REDUCTION ARE STORED IN THE EIS3
+C REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. EIS3
+C EIS3
+C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS EIS3
+C INTERCHANGED IN THE REDUCTION. EIS3
+C ONLY ELEMENTS LOW THROUGH IGH ARE USED. EIS3
+C EIS3
+C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, EIS3
+C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY EIS3
+C EIS3
+C THIS VERSION DATED AUGUST 1983. EIS3
+C EIS3
+C ------------------------------------------------------------------EIS3
+C EIS3
+ LA = IGH - 1 EIS3
+ KP1 = LOW + 1 EIS3
+ IF (LA .LT. KP1) GO TO 200 EIS3
+C EIS3
+ DO 180 M = KP1, LA EIS3
+ MM1 = M - 1 EIS3
+ X = 0.0D0 EIS3
+ I = M EIS3
+C EIS3
+ DO 100 J = M, IGH EIS3
+ IF (DABS(A(J,MM1)) .LE. DABS(X)) GO TO 100 EIS3
+ X = A(J,MM1) EIS3
+ I = J EIS3
+ 100 CONTINUE EIS3
+C EIS3
+ INT(M) = I EIS3
+ IF (I .EQ. M) GO TO 130 EIS3
+C .......... INTERCHANGE ROWS AND COLUMNS OF A .......... EIS3
+ DO 110 J = MM1, N EIS3
+ Y = A(I,J) EIS3
+ A(I,J) = A(M,J) EIS3
+ A(M,J) = Y EIS3
+ 110 CONTINUE EIS3
+C EIS3
+ DO 120 J = 1, IGH EIS3
+ Y = A(J,I) EIS3
+ A(J,I) = A(J,M) EIS3
+ A(J,M) = Y EIS3
+ 120 CONTINUE EIS3
+C .......... END INTERCHANGE .......... EIS3
+ 130 IF (X .EQ. 0.0D0) GO TO 180 EIS3
+ MP1 = M + 1 EIS3
+C EIS3
+ DO 160 I = MP1, IGH EIS3
+ Y = A(I,MM1) EIS3
+ IF (Y .EQ. 0.0D0) GO TO 160 EIS3
+ Y = Y / X EIS3
+ A(I,MM1) = Y EIS3
+C EIS3
+ DO 140 J = M, N EIS3
+ 140 A(I,J) = A(I,J) - Y * A(M,J) EIS3
+C EIS3
+ DO 150 J = 1, IGH EIS3
+ 150 A(J,M) = A(J,M) + Y * A(J,I) EIS3
+C EIS3
+ 160 CONTINUE EIS3
+C EIS3
+ 180 CONTINUE EIS3
+C EIS3
+ 200 RETURN EIS3
+ END EIS3
Added: grass-addons/m.eigensystem/eltran.f
===================================================================
--- grass-addons/m.eigensystem/eltran.f (rev 0)
+++ grass-addons/m.eigensystem/eltran.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,78 @@
+ SUBROUTINE ELTRAN(NM,N,LOW,IGH,A,INT,Z) EIS4
+C EIS4
+ INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1 EIS4
+ DOUBLE PRECISION A(NM,IGH),Z(NM,N) EIS4
+ INTEGER INT(IGH) EIS4
+C EIS4
+C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS, EIS4
+C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. EIS4
+C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). EIS4
+C EIS4
+C THIS SUBROUTINE ACCUMULATES THE STABILIZED ELEMENTARY EIS4
+C SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A EIS4
+C REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY ELMHES. EIS4
+C EIS4
+C ON INPUT EIS4
+C EIS4
+C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL EIS4
+C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM EIS4
+C DIMENSION STATEMENT. EIS4
+C EIS4
+C N IS THE ORDER OF THE MATRIX. EIS4
+C EIS4
+C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING EIS4
+C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, EIS4
+C SET LOW=1, IGH=N. EIS4
+C EIS4
+C A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE EIS4
+C REDUCTION BY ELMHES IN ITS LOWER TRIANGLE EIS4
+C BELOW THE SUBDIAGONAL. EIS4
+C EIS4
+C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS EIS4
+C INTERCHANGED IN THE REDUCTION BY ELMHES. EIS4
+C ONLY ELEMENTS LOW THROUGH IGH ARE USED. EIS4
+C EIS4
+C ON OUTPUT EIS4
+C EIS4
+C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE EIS4
+C REDUCTION BY ELMHES. EIS4
+C EIS4
+C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, EIS4
+C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY EIS4
+C EIS4
+C THIS VERSION DATED AUGUST 1983. EIS4
+C EIS4
+C ------------------------------------------------------------------EIS4
+C EIS4
+C .......... INITIALIZE Z TO IDENTITY MATRIX .......... EIS4
+ DO 80 J = 1, N EIS4
+C EIS4
+ DO 60 I = 1, N EIS4
+ 60 Z(I,J) = 0.0D0 EIS4
+C EIS4
+ Z(J,J) = 1.0D0 EIS4
+ 80 CONTINUE EIS4
+C EIS4
+ KL = IGH - LOW - 1 EIS4
+ IF (KL .LT. 1) GO TO 200 EIS4
+C .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... EIS4
+ DO 140 MM = 1, KL EIS4
+ MP = IGH - MM EIS4
+ MP1 = MP + 1 EIS4
+C EIS4
+ DO 100 I = MP1, IGH EIS4
+ 100 Z(I,MP) = A(I,MP-1) EIS4
+C EIS4
+ I = INT(MP) EIS4
+ IF (I .EQ. MP) GO TO 140 EIS4
+C EIS4
+ DO 130 J = MP, IGH EIS4
+ Z(MP,J) = Z(I,J) EIS4
+ Z(I,J) = 0.0D0 EIS4
+ 130 CONTINUE EIS4
+C EIS4
+ Z(I,MP) = 1.0D0 EIS4
+ 140 CONTINUE EIS4
+C EIS4
+ 200 RETURN EIS4
+ END EIS4
Added: grass-addons/m.eigensystem/hqr.f
===================================================================
--- grass-addons/m.eigensystem/hqr.f (rev 0)
+++ grass-addons/m.eigensystem/hqr.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,234 @@
+ SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR) EIS4
+C EIS4
+ INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR EIS4
+ DOUBLE PRECISION H(NM,N),WR(N),WI(N) EIS4
+ DOUBLE PRECISION P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2 EIS4
+ LOGICAL NOTLAS EIS4
+C EIS4
+C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, EIS4
+C NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. EIS4
+C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). EIS4
+C EIS4
+C THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL EIS4
+C UPPER HESSENBERG MATRIX BY THE QR METHOD. EIS4
+C EIS4
+C ON INPUT EIS4
+C EIS4
+C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL EIS4
+C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM EIS4
+C DIMENSION STATEMENT. EIS4
+C EIS4
+C N IS THE ORDER OF THE MATRIX. EIS4
+C EIS4
+C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING EIS4
+C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, EIS4
+C SET LOW=1, IGH=N. EIS4
+C EIS4
+C H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT EIS4
+C THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG EIS4
+C FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED EIS4
+C IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. EIS4
+C EIS4
+C ON OUTPUT EIS4
+C EIS4
+C H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED EIS4
+C BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND EIS4
+C BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. EIS4
+C EIS4
+C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, EIS4
+C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES EIS4
+C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS EIS4
+C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE EIS4
+C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN EIS4
+C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT EIS4
+C FOR INDICES IERR+1,...,N. EIS4
+C EIS4
+C IERR IS SET TO EIS4
+C ZERO FOR NORMAL RETURN, EIS4
+C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED EIS4
+C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. EIS4
+C EIS4
+C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, EIS4
+C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY EIS4
+C EIS4
+C THIS VERSION DATED AUGUST 1983. EIS4
+C EIS4
+C ------------------------------------------------------------------EIS4
+C EIS4
+ IERR = 0 EIS4
+ NORM = 0.0D0 EIS4
+ K = 1 EIS4
+C .......... STORE ROOTS ISOLATED BY BALANC EIS4
+C AND COMPUTE MATRIX NORM .......... EIS4
+ DO 50 I = 1, N EIS4
+C EIS4
+ DO 40 J = K, N EIS4
+ 40 NORM = NORM + DABS(H(I,J)) EIS4
+C EIS4
+ K = I EIS4
+ IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 EIS4
+ WR(I) = H(I,I) EIS4
+ WI(I) = 0.0D0 EIS4
+ 50 CONTINUE EIS4
+C EIS4
+ EN = IGH EIS4
+ T = 0.0D0 EIS4
+ ITN = 30*N EIS4
+C .......... SEARCH FOR NEXT EIGENVALUES .......... EIS4
+ 60 IF (EN .LT. LOW) GO TO 1001 EIS4
+ ITS = 0 EIS4
+ NA = EN - 1 EIS4
+ ENM2 = NA - 1 EIS4
+C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT EIS4
+C FOR L=EN STEP -1 UNTIL LOW DO -- .......... EIS4
+ 70 DO 80 LL = LOW, EN EIS4
+ L = EN + LOW - LL EIS4
+ IF (L .EQ. LOW) GO TO 100 EIS4
+ S = DABS(H(L-1,L-1)) + DABS(H(L,L)) EIS4
+ IF (S .EQ. 0.0D0) S = NORM EIS4
+ TST1 = S EIS4
+ TST2 = TST1 + DABS(H(L,L-1)) EIS4
+ IF (TST2 .EQ. TST1) GO TO 100 EIS4
+ 80 CONTINUE EIS4
+C .......... FORM SHIFT .......... EIS4
+ 100 X = H(EN,EN) EIS4
+ IF (L .EQ. EN) GO TO 270 EIS4
+ Y = H(NA,NA) EIS4
+ W = H(EN,NA) * H(NA,EN) EIS4
+ IF (L .EQ. NA) GO TO 280 EIS4
+ IF (ITN .EQ. 0) GO TO 1000 EIS4
+ IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 EIS4
+C .......... FORM EXCEPTIONAL SHIFT .......... EIS4
+ T = T + X EIS4
+C EIS4
+ DO 120 I = LOW, EN EIS4
+ 120 H(I,I) = H(I,I) - X EIS4
+C EIS4
+ S = DABS(H(EN,NA)) + DABS(H(NA,ENM2)) EIS4
+ X = 0.75D0 * S EIS4
+ Y = X EIS4
+ W = -0.4375D0 * S * S EIS4
+ 130 ITS = ITS + 1 EIS4
+ ITN = ITN - 1 EIS4
+C .......... LOOK FOR TWO CONSECUTIVE SMALL EIS4
+C SUB-DIAGONAL ELEMENTS. EIS4
+C FOR M=EN-2 STEP -1 UNTIL L DO -- .......... EIS4
+ DO 140 MM = L, ENM2 EIS4
+ M = ENM2 + L - MM EIS4
+ ZZ = H(M,M) EIS4
+ R = X - ZZ EIS4
+ S = Y - ZZ EIS4
+ P = (R * S - W) / H(M+1,M) + H(M,M+1) EIS4
+ Q = H(M+1,M+1) - ZZ - R - S EIS4
+ R = H(M+2,M+1) EIS4
+ S = DABS(P) + DABS(Q) + DABS(R) EIS4
+ P = P / S EIS4
+ Q = Q / S EIS4
+ R = R / S EIS4
+ IF (M .EQ. L) GO TO 150 EIS4
+ TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))EIS4
+ TST2 = TST1 + DABS(H(M,M-1))*(DABS(Q) + DABS(R)) EIS4
+ IF (TST2 .EQ. TST1) GO TO 150 EIS4
+ 140 CONTINUE EIS4
+C EIS4
+ 150 MP2 = M + 2 EIS4
+C EIS4
+ DO 160 I = MP2, EN EIS4
+ H(I,I-2) = 0.0D0 EIS4
+ IF (I .EQ. MP2) GO TO 160 EIS4
+ H(I,I-3) = 0.0D0 EIS4
+ 160 CONTINUE EIS4
+C .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND EIS4
+C COLUMNS M TO EN .......... EIS4
+ DO 260 K = M, NA EIS4
+ NOTLAS = K .NE. NA EIS4
+ IF (K .EQ. M) GO TO 170 EIS4
+ P = H(K,K-1) EIS4
+ Q = H(K+1,K-1) EIS4
+ R = 0.0D0 EIS4
+ IF (NOTLAS) R = H(K+2,K-1) EIS4
+ X = DABS(P) + DABS(Q) + DABS(R) EIS4
+ IF (X .EQ. 0.0D0) GO TO 260 EIS4
+ P = P / X EIS4
+ Q = Q / X EIS4
+ R = R / X EIS4
+ 170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P) EIS4
+ IF (K .EQ. M) GO TO 180 EIS4
+ H(K,K-1) = -S * X EIS4
+ GO TO 190 EIS4
+ 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) EIS4
+ 190 P = P + S EIS4
+ X = P / S EIS4
+ Y = Q / S EIS4
+ ZZ = R / S EIS4
+ Q = Q / P EIS4
+ R = R / P EIS4
+ IF (NOTLAS) GO TO 225 EIS4
+C .......... ROW MODIFICATION .......... EIS4
+ DO 200 J = K, N EIS4
+ P = H(K,J) + Q * H(K+1,J) EIS4
+ H(K,J) = H(K,J) - P * X EIS4
+ H(K+1,J) = H(K+1,J) - P * Y EIS4
+ 200 CONTINUE EIS4
+C EIS4
+ J = MIN0(EN,K+3) EIS4
+C .......... COLUMN MODIFICATION .......... EIS4
+ DO 210 I = 1, J EIS4
+ P = X * H(I,K) + Y * H(I,K+1) EIS4
+ H(I,K) = H(I,K) - P EIS4
+ H(I,K+1) = H(I,K+1) - P * Q EIS4
+ 210 CONTINUE EIS4
+ GO TO 255 EIS4
+ 225 CONTINUE EIS4
+C .......... ROW MODIFICATION .......... EIS4
+ DO 230 J = K, N EIS4
+ P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J) EIS4
+ H(K,J) = H(K,J) - P * X EIS4
+ H(K+1,J) = H(K+1,J) - P * Y EIS4
+ H(K+2,J) = H(K+2,J) - P * ZZ EIS4
+ 230 CONTINUE EIS4
+C EIS4
+ J = MIN0(EN,K+3) EIS4
+C .......... COLUMN MODIFICATION .......... EIS4
+ DO 240 I = 1, J EIS4
+ P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2) EIS4
+ H(I,K) = H(I,K) - P EIS4
+ H(I,K+1) = H(I,K+1) - P * Q EIS4
+ H(I,K+2) = H(I,K+2) - P * R EIS4
+ 240 CONTINUE EIS4
+ 255 CONTINUE EIS4
+C EIS4
+ 260 CONTINUE EIS4
+C EIS4
+ GO TO 70 EIS4
+C .......... ONE ROOT FOUND .......... EIS4
+ 270 WR(EN) = X + T EIS4
+ WI(EN) = 0.0D0 EIS4
+ EN = NA EIS4
+ GO TO 60 EIS4
+C .......... TWO ROOTS FOUND .......... EIS4
+ 280 P = (Y - X) / 2.0D0 EIS4
+ Q = P * P + W EIS4
+ ZZ = DSQRT(DABS(Q)) EIS4
+ X = X + T EIS4
+ IF (Q .LT. 0.0D0) GO TO 320 EIS4
+C .......... REAL PAIR .......... EIS4
+ ZZ = P + DSIGN(ZZ,P) EIS4
+ WR(NA) = X + ZZ EIS4
+ WR(EN) = WR(NA) EIS4
+ IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ EIS4
+ WI(NA) = 0.0D0 EIS4
+ WI(EN) = 0.0D0 EIS4
+ GO TO 330 EIS4
+C .......... COMPLEX PAIR .......... EIS4
+ 320 WR(NA) = X + P EIS4
+ WR(EN) = X + P EIS4
+ WI(NA) = ZZ EIS4
+ WI(EN) = -ZZ EIS4
+ 330 EN = ENM2 EIS4
+ GO TO 60 EIS4
+C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT EIS4
+C CONVERGED AFTER 30*N ITERATIONS .......... EIS4
+ 1000 IERR = EN EIS4
+ 1001 RETURN EIS4
+ END EIS4
Added: grass-addons/m.eigensystem/hqr2.f
===================================================================
--- grass-addons/m.eigensystem/hqr2.f (rev 0)
+++ grass-addons/m.eigensystem/hqr2.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,449 @@
+ SUBROUTINE HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR) EIS4
+C EIS4
+ INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN, EIS4
+ X IGH,ITN,ITS,LOW,MP2,ENM2,IERR EIS4
+ DOUBLE PRECISION H(NM,N),WR(N),WI(N),Z(NM,N) EIS4
+ DOUBLE PRECISION P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2 EIS4
+ LOGICAL NOTLAS EIS4
+C EIS4
+C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, EIS4
+C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. EIS4
+C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). EIS4
+C EIS4
+C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS EIS4
+C OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD. THE EIS4
+C EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND EIS4
+C IF ELMHES AND ELTRAN OR ORTHES AND ORTRAN HAVE EIS4
+C BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM EIS4
+C AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS. EIS4
+C EIS4
+C ON INPUT EIS4
+C EIS4
+C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL EIS4
+C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM EIS4
+C DIMENSION STATEMENT. EIS4
+C EIS4
+C N IS THE ORDER OF THE MATRIX. EIS4
+C EIS4
+C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING EIS4
+C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, EIS4
+C SET LOW=1, IGH=N. EIS4
+C EIS4
+C H CONTAINS THE UPPER HESSENBERG MATRIX. EIS4
+C EIS4
+C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ELTRAN EIS4
+C AFTER THE REDUCTION BY ELMHES, OR BY ORTRAN AFTER THE EIS4
+C REDUCTION BY ORTHES, IF PERFORMED. IF THE EIGENVECTORS EIS4
+C OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE EIS4
+C IDENTITY MATRIX. EIS4
+C EIS4
+C ON OUTPUT EIS4
+C EIS4
+C H HAS BEEN DESTROYED. EIS4
+C EIS4
+C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, EIS4
+C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES EIS4
+C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS EIS4
+C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE EIS4
+C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN EIS4
+C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT EIS4
+C FOR INDICES IERR+1,...,N. EIS4
+C EIS4
+C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. EIS4
+C IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z EIS4
+C CONTAINS ITS EIGENVECTOR. IF THE I-TH EIGENVALUE IS COMPLEX EIS4
+C WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH EIS4
+C COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS EIS4
+C EIGENVECTOR. THE EIGENVECTORS ARE UNNORMALIZED. IF AN EIS4
+C ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND. EIS4
+C EIS4
+C IERR IS SET TO EIS4
+C ZERO FOR NORMAL RETURN, EIS4
+C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED EIS4
+C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. EIS4
+C EIS4
+C CALLS CDIV FOR COMPLEX DIVISION. EIS4
+C EIS4
+C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, EIS4
+C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY EIS4
+C EIS4
+C THIS VERSION DATED AUGUST 1983. EIS4
+C EIS4
+C ------------------------------------------------------------------EIS4
+C EIS4
+ IERR = 0 EIS4
+ NORM = 0.0D0 EIS4
+ K = 1 EIS4
+C .......... STORE ROOTS ISOLATED BY BALANC EIS4
+C AND COMPUTE MATRIX NORM .......... EIS4
+ DO 50 I = 1, N EIS4
+C EIS4
+ DO 40 J = K, N EIS4
+ 40 NORM = NORM + DABS(H(I,J)) EIS4
+C EIS4
+ K = I EIS4
+ IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 EIS4
+ WR(I) = H(I,I) EIS4
+ WI(I) = 0.0D0 EIS4
+ 50 CONTINUE EIS4
+C EIS4
+ EN = IGH EIS4
+ T = 0.0D0 EIS4
+ ITN = 30*N EIS4
+C .......... SEARCH FOR NEXT EIGENVALUES .......... EIS4
+ 60 IF (EN .LT. LOW) GO TO 340 EIS4
+ ITS = 0 EIS4
+ NA = EN - 1 EIS4
+ ENM2 = NA - 1 EIS4
+C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT EIS4
+C FOR L=EN STEP -1 UNTIL LOW DO -- .......... EIS4
+ 70 DO 80 LL = LOW, EN EIS4
+ L = EN + LOW - LL EIS4
+ IF (L .EQ. LOW) GO TO 100 EIS4
+ S = DABS(H(L-1,L-1)) + DABS(H(L,L)) EIS4
+ IF (S .EQ. 0.0D0) S = NORM EIS4
+ TST1 = S EIS4
+ TST2 = TST1 + DABS(H(L,L-1)) EIS4
+ IF (TST2 .EQ. TST1) GO TO 100 EIS4
+ 80 CONTINUE EIS4
+C .......... FORM SHIFT .......... EIS4
+ 100 X = H(EN,EN) EIS4
+ IF (L .EQ. EN) GO TO 270 EIS4
+ Y = H(NA,NA) EIS4
+ W = H(EN,NA) * H(NA,EN) EIS4
+ IF (L .EQ. NA) GO TO 280 EIS4
+ IF (ITN .EQ. 0) GO TO 1000 EIS4
+ IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 EIS4
+C .......... FORM EXCEPTIONAL SHIFT .......... EIS4
+ T = T + X EIS4
+C EIS4
+ DO 120 I = LOW, EN EIS4
+ 120 H(I,I) = H(I,I) - X EIS4
+C EIS4
+ S = DABS(H(EN,NA)) + DABS(H(NA,ENM2)) EIS4
+ X = 0.75D0 * S EIS4
+ Y = X EIS4
+ W = -0.4375D0 * S * S EIS4
+ 130 ITS = ITS + 1 EIS4
+ ITN = ITN - 1 EIS4
+C .......... LOOK FOR TWO CONSECUTIVE SMALL EIS4
+C SUB-DIAGONAL ELEMENTS. EIS4
+C FOR M=EN-2 STEP -1 UNTIL L DO -- .......... EIS4
+ DO 140 MM = L, ENM2 EIS4
+ M = ENM2 + L - MM EIS4
+ ZZ = H(M,M) EIS4
+ R = X - ZZ EIS4
+ S = Y - ZZ EIS4
+ P = (R * S - W) / H(M+1,M) + H(M,M+1) EIS4
+ Q = H(M+1,M+1) - ZZ - R - S EIS4
+ R = H(M+2,M+1) EIS4
+ S = DABS(P) + DABS(Q) + DABS(R) EIS4
+ P = P / S EIS4
+ Q = Q / S EIS4
+ R = R / S EIS4
+ IF (M .EQ. L) GO TO 150 EIS4
+ TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))EIS4
+ TST2 = TST1 + DABS(H(M,M-1))*(DABS(Q) + DABS(R)) EIS4
+ IF (TST2 .EQ. TST1) GO TO 150 EIS4
+ 140 CONTINUE EIS4
+C EIS4
+ 150 MP2 = M + 2 EIS4
+C EIS4
+ DO 160 I = MP2, EN EIS4
+ H(I,I-2) = 0.0D0 EIS4
+ IF (I .EQ. MP2) GO TO 160 EIS4
+ H(I,I-3) = 0.0D0 EIS4
+ 160 CONTINUE EIS4
+C .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND EIS4
+C COLUMNS M TO EN .......... EIS4
+ DO 260 K = M, NA EIS4
+ NOTLAS = K .NE. NA EIS4
+ IF (K .EQ. M) GO TO 170 EIS4
+ P = H(K,K-1) EIS4
+ Q = H(K+1,K-1) EIS4
+ R = 0.0D0 EIS4
+ IF (NOTLAS) R = H(K+2,K-1) EIS4
+ X = DABS(P) + DABS(Q) + DABS(R) EIS4
+ IF (X .EQ. 0.0D0) GO TO 260 EIS4
+ P = P / X EIS4
+ Q = Q / X EIS4
+ R = R / X EIS4
+ 170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P) EIS4
+ IF (K .EQ. M) GO TO 180 EIS4
+ H(K,K-1) = -S * X EIS4
+ GO TO 190 EIS4
+ 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) EIS4
+ 190 P = P + S EIS4
+ X = P / S EIS4
+ Y = Q / S EIS4
+ ZZ = R / S EIS4
+ Q = Q / P EIS4
+ R = R / P EIS4
+ IF (NOTLAS) GO TO 225 EIS4
+C .......... ROW MODIFICATION .......... EIS4
+ DO 200 J = K, N EIS4
+ P = H(K,J) + Q * H(K+1,J) EIS4
+ H(K,J) = H(K,J) - P * X EIS4
+ H(K+1,J) = H(K+1,J) - P * Y EIS4
+ 200 CONTINUE EIS4
+C EIS4
+ J = MIN0(EN,K+3) EIS4
+C .......... COLUMN MODIFICATION .......... EIS4
+ DO 210 I = 1, J EIS4
+ P = X * H(I,K) + Y * H(I,K+1) EIS4
+ H(I,K) = H(I,K) - P EIS4
+ H(I,K+1) = H(I,K+1) - P * Q EIS4
+ 210 CONTINUE EIS4
+C .......... ACCUMULATE TRANSFORMATIONS .......... EIS4
+ DO 220 I = LOW, IGH EIS4
+ P = X * Z(I,K) + Y * Z(I,K+1) EIS4
+ Z(I,K) = Z(I,K) - P EIS4
+ Z(I,K+1) = Z(I,K+1) - P * Q EIS4
+ 220 CONTINUE EIS4
+ GO TO 255 EIS4
+ 225 CONTINUE EIS4
+C .......... ROW MODIFICATION .......... EIS4
+ DO 230 J = K, N EIS4
+ P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J) EIS4
+ H(K,J) = H(K,J) - P * X EIS4
+ H(K+1,J) = H(K+1,J) - P * Y EIS4
+ H(K+2,J) = H(K+2,J) - P * ZZ EIS4
+ 230 CONTINUE EIS4
+C EIS4
+ J = MIN0(EN,K+3) EIS4
+C .......... COLUMN MODIFICATION .......... EIS4
+ DO 240 I = 1, J EIS4
+ P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2) EIS4
+ H(I,K) = H(I,K) - P EIS4
+ H(I,K+1) = H(I,K+1) - P * Q EIS4
+ H(I,K+2) = H(I,K+2) - P * R EIS4
+ 240 CONTINUE EIS4
+C .......... ACCUMULATE TRANSFORMATIONS .......... EIS4
+ DO 250 I = LOW, IGH EIS4
+ P = X * Z(I,K) + Y * Z(I,K+1) + ZZ * Z(I,K+2) EIS4
+ Z(I,K) = Z(I,K) - P EIS4
+ Z(I,K+1) = Z(I,K+1) - P * Q EIS4
+ Z(I,K+2) = Z(I,K+2) - P * R EIS4
+ 250 CONTINUE EIS4
+ 255 CONTINUE EIS4
+C EIS4
+ 260 CONTINUE EIS4
+C EIS4
+ GO TO 70 EIS4
+C .......... ONE ROOT FOUND .......... EIS4
+ 270 H(EN,EN) = X + T EIS4
+ WR(EN) = H(EN,EN) EIS4
+ WI(EN) = 0.0D0 EIS4
+ EN = NA EIS4
+ GO TO 60 EIS4
+C .......... TWO ROOTS FOUND .......... EIS4
+ 280 P = (Y - X) / 2.0D0 EIS4
+ Q = P * P + W EIS4
+ ZZ = DSQRT(DABS(Q)) EIS4
+ H(EN,EN) = X + T EIS4
+ X = H(EN,EN) EIS4
+ H(NA,NA) = Y + T EIS4
+ IF (Q .LT. 0.0D0) GO TO 320 EIS4
+C .......... REAL PAIR .......... EIS4
+ ZZ = P + DSIGN(ZZ,P) EIS4
+ WR(NA) = X + ZZ EIS4
+ WR(EN) = WR(NA) EIS4
+ IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ EIS4
+ WI(NA) = 0.0D0 EIS4
+ WI(EN) = 0.0D0 EIS4
+ X = H(EN,NA) EIS4
+ S = DABS(X) + DABS(ZZ) EIS4
+ P = X / S EIS4
+ Q = ZZ / S EIS4
+ R = DSQRT(P*P+Q*Q) EIS4
+ P = P / R EIS4
+ Q = Q / R EIS4
+C .......... ROW MODIFICATION .......... EIS4
+ DO 290 J = NA, N EIS4
+ ZZ = H(NA,J) EIS4
+ H(NA,J) = Q * ZZ + P * H(EN,J) EIS4
+ H(EN,J) = Q * H(EN,J) - P * ZZ EIS4
+ 290 CONTINUE EIS4
+C .......... COLUMN MODIFICATION .......... EIS4
+ DO 300 I = 1, EN EIS4
+ ZZ = H(I,NA) EIS4
+ H(I,NA) = Q * ZZ + P * H(I,EN) EIS4
+ H(I,EN) = Q * H(I,EN) - P * ZZ EIS4
+ 300 CONTINUE EIS4
+C .......... ACCUMULATE TRANSFORMATIONS .......... EIS4
+ DO 310 I = LOW, IGH EIS4
+ ZZ = Z(I,NA) EIS4
+ Z(I,NA) = Q * ZZ + P * Z(I,EN) EIS4
+ Z(I,EN) = Q * Z(I,EN) - P * ZZ EIS4
+ 310 CONTINUE EIS4
+C EIS4
+ GO TO 330 EIS4
+C .......... COMPLEX PAIR .......... EIS4
+ 320 WR(NA) = X + P EIS4
+ WR(EN) = X + P EIS4
+ WI(NA) = ZZ EIS4
+ WI(EN) = -ZZ EIS4
+ 330 EN = ENM2 EIS4
+ GO TO 60 EIS4
+C .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND EIS4
+C VECTORS OF UPPER TRIANGULAR FORM .......... EIS4
+ 340 IF (NORM .EQ. 0.0D0) GO TO 1001 EIS4
+C .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... EIS4
+ DO 800 NN = 1, N EIS4
+ EN = N + 1 - NN EIS4
+ P = WR(EN) EIS4
+ Q = WI(EN) EIS4
+ NA = EN - 1 EIS4
+ IF (Q) 710, 600, 800 EIS4
+C .......... REAL VECTOR .......... EIS4
+ 600 M = EN EIS4
+ H(EN,EN) = 1.0D0 EIS4
+ IF (NA .EQ. 0) GO TO 800 EIS4
+C .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... EIS4
+ DO 700 II = 1, NA EIS4
+ I = EN - II EIS4
+ W = H(I,I) - P EIS4
+ R = 0.0D0 EIS4
+C EIS4
+ DO 610 J = M, EN EIS4
+ 610 R = R + H(I,J) * H(J,EN) EIS4
+C EIS4
+ IF (WI(I) .GE. 0.0D0) GO TO 630 EIS4
+ ZZ = W EIS4
+ S = R EIS4
+ GO TO 700 EIS4
+ 630 M = I EIS4
+ IF (WI(I) .NE. 0.0D0) GO TO 640 EIS4
+ T = W EIS4
+ IF (T .NE. 0.0D0) GO TO 635 EIS4
+ TST1 = NORM EIS4
+ T = TST1 EIS4
+ 632 T = 0.01D0 * T EIS4
+ TST2 = NORM + T EIS4
+ IF (TST2 .GT. TST1) GO TO 632 EIS4
+ 635 H(I,EN) = -R / T EIS4
+ GO TO 680 EIS4
+C .......... SOLVE REAL EQUATIONS .......... EIS4
+ 640 X = H(I,I+1) EIS4
+ Y = H(I+1,I) EIS4
+ Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) EIS4
+ T = (X * S - ZZ * R) / Q EIS4
+ H(I,EN) = T EIS4
+ IF (DABS(X) .LE. DABS(ZZ)) GO TO 650 EIS4
+ H(I+1,EN) = (-R - W * T) / X EIS4
+ GO TO 680 EIS4
+ 650 H(I+1,EN) = (-S - Y * T) / ZZ EIS4
+C EIS4
+C .......... OVERFLOW CONTROL .......... EIS4
+ 680 T = DABS(H(I,EN)) EIS4
+ IF (T .EQ. 0.0D0) GO TO 700 EIS4
+ TST1 = T EIS4
+ TST2 = TST1 + 1.0D0/TST1 EIS4
+ IF (TST2 .GT. TST1) GO TO 700 EIS4
+ DO 690 J = I, EN EIS4
+ H(J,EN) = H(J,EN)/T EIS4
+ 690 CONTINUE EIS4
+C EIS4
+ 700 CONTINUE EIS4
+C .......... END REAL VECTOR .......... EIS4
+ GO TO 800 EIS4
+C .......... COMPLEX VECTOR .......... EIS4
+ 710 M = NA EIS4
+C .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT EIS4
+C EIGENVECTOR MATRIX IS TRIANGULAR .......... EIS4
+ IF (DABS(H(EN,NA)) .LE. DABS(H(NA,EN))) GO TO 720 EIS4
+ H(NA,NA) = Q / H(EN,NA) EIS4
+ H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA) EIS4
+ GO TO 730 EIS4
+ 720 CALL CDIV(0.0D0,-H(NA,EN),H(NA,NA)-P,Q,H(NA,NA),H(NA,EN)) EIS4
+ 730 H(EN,NA) = 0.0D0 EIS4
+ H(EN,EN) = 1.0D0 EIS4
+ ENM2 = NA - 1 EIS4
+ IF (ENM2 .EQ. 0) GO TO 800 EIS4
+C .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... EIS4
+ DO 795 II = 1, ENM2 EIS4
+ I = NA - II EIS4
+ W = H(I,I) - P EIS4
+ RA = 0.0D0 EIS4
+ SA = 0.0D0 EIS4
+C EIS4
+ DO 760 J = M, EN EIS4
+ RA = RA + H(I,J) * H(J,NA) EIS4
+ SA = SA + H(I,J) * H(J,EN) EIS4
+ 760 CONTINUE EIS4
+C EIS4
+ IF (WI(I) .GE. 0.0D0) GO TO 770 EIS4
+ ZZ = W EIS4
+ R = RA EIS4
+ S = SA EIS4
+ GO TO 795 EIS4
+ 770 M = I EIS4
+ IF (WI(I) .NE. 0.0D0) GO TO 780 EIS4
+ CALL CDIV(-RA,-SA,W,Q,H(I,NA),H(I,EN)) EIS4
+ GO TO 790 EIS4
+C .......... SOLVE COMPLEX EQUATIONS .......... EIS4
+ 780 X = H(I,I+1) EIS4
+ Y = H(I+1,I) EIS4
+ VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q EIS4
+ VI = (WR(I) - P) * 2.0D0 * Q EIS4
+ IF (VR .NE. 0.0D0 .OR. VI .NE. 0.0D0) GO TO 784 EIS4
+ TST1 = NORM * (DABS(W) + DABS(Q) + DABS(X) EIS4
+ X + DABS(Y) + DABS(ZZ)) EIS4
+ VR = TST1 EIS4
+ 783 VR = 0.01D0 * VR EIS4
+ TST2 = TST1 + VR EIS4
+ IF (TST2 .GT. TST1) GO TO 783 EIS4
+ 784 CALL CDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI, EIS4
+ X H(I,NA),H(I,EN)) EIS4
+ IF (DABS(X) .LE. DABS(ZZ) + DABS(Q)) GO TO 785 EIS4
+ H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X EIS4
+ H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X EIS4
+ GO TO 790 EIS4
+ 785 CALL CDIV(-R-Y*H(I,NA),-S-Y*H(I,EN),ZZ,Q, EIS4
+ X H(I+1,NA),H(I+1,EN)) EIS4
+C EIS4
+C .......... OVERFLOW CONTROL .......... EIS4
+ 790 T = DMAX1(DABS(H(I,NA)), DABS(H(I,EN))) EIS4
+ IF (T .EQ. 0.0D0) GO TO 795 EIS4
+ TST1 = T EIS4
+ TST2 = TST1 + 1.0D0/TST1 EIS4
+ IF (TST2 .GT. TST1) GO TO 795 EIS4
+ DO 792 J = I, EN EIS4
+ H(J,NA) = H(J,NA)/T EIS4
+ H(J,EN) = H(J,EN)/T EIS4
+ 792 CONTINUE EIS4
+C EIS4
+ 795 CONTINUE EIS4
+C .......... END COMPLEX VECTOR .......... EIS4
+ 800 CONTINUE EIS4
+C .......... END BACK SUBSTITUTION. EIS4
+C VECTORS OF ISOLATED ROOTS .......... EIS4
+ DO 840 I = 1, N EIS4
+ IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840 EIS4
+C EIS4
+ DO 820 J = I, N EIS4
+ 820 Z(I,J) = H(I,J) EIS4
+C EIS4
+ 840 CONTINUE EIS4
+C .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE EIS4
+C VECTORS OF ORIGINAL FULL MATRIX. EIS4
+C FOR J=N STEP -1 UNTIL LOW DO -- .......... EIS4
+ DO 880 JJ = LOW, N EIS4
+ J = N + LOW - JJ EIS4
+ M = MIN0(J,IGH) EIS4
+C EIS4
+ DO 880 I = LOW, IGH EIS4
+ ZZ = 0.0D0 EIS4
+C EIS4
+ DO 860 K = LOW, M EIS4
+ 860 ZZ = ZZ + Z(I,K) * H(K,J) EIS4
+C EIS4
+ Z(I,J) = ZZ EIS4
+ 880 CONTINUE EIS4
+C EIS4
+ GO TO 1001 EIS4
+C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT EIS4
+C CONVERGED AFTER 30*N ITERATIONS .......... EIS4
+ 1000 IERR = EN EIS4
+ 1001 RETURN EIS4
+ END EIS4
Added: grass-addons/m.eigensystem/m.eigensystem.html
===================================================================
--- grass-addons/m.eigensystem/m.eigensystem.html (rev 0)
+++ grass-addons/m.eigensystem/m.eigensystem.html 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,180 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS: m.eigensystem</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<H2>NAME</H2>
+
+<EM><b>m.eigensystem</b></EM> - Computes eigen values and eigen vectors
+for square matricies.
+
+<BR>
+
+<EM>(GRASS Data Import/Processing Program)</EM>
+
+<H2>SYNOPSIS</H2>
+
+<B>m.eigensystem < </B><EM>inputfile </EM>
+
+<H2>DESCRIPTION</H2>
+
+<EM>m.eigensystem</EM>
+determines the eigen values and eigen vectors for square matricies.
+The
+<EM>inputfile</EM>
+must have the following format:
+the first line contains an integer K
+which is the number of rows and columns in the matrix;
+the remainder of the file is the matrix, i.e.,
+K lines, each containing K real numbers.
+For example:
+<P>
+<PRE>
+ 3
+ 462.876649 480.411218 281.758307
+ 480.411218 513.015646 278.914813
+ 281.758307 278.914813 336.326645
+
+</PRE>
+<P>
+
+The output will be K groups of lines; each group will have the format:
+<PRE>
+ E real part imaginary part relative importance
+ V real part imaginary part
+ ... K lines ...
+ N real part imaginary part
+ ... K lines ...
+ W real part imaginary part
+ ... K lines ...
+
+</PRE>
+
+The
+<EM>E</EM>
+line is the eigen value.
+The
+<EM>V</EM>
+lines are the eigen vector associated with E.
+The
+<EM>N</EM>
+lines are the V vector normalized to have a magnitude of 1.
+The
+<EM>W</EM>
+lines are the N vector multiplied by the square root of the
+magnitude of the eigen value (E).
+
+
+<P>
+
+For the example input matrix, the output would be:
+<PRE>
+ E 1159.7452017844 0.0000000000 88.38
+ V 0.6910021591 0.0000000000
+ V 0.7205280412 0.0000000000
+ V 0.4805108400 0.0000000000
+ N 0.6236808478 0.0000000000
+ N 0.6503301526 0.0000000000
+ N 0.4336967751 0.0000000000
+ W 21.2394712045 0.0000000000
+ W 22.1470141296 0.0000000000
+ W 14.7695575384 0.0000000000
+
+ E 5.9705414972 0.0000000000 0.45
+ V 0.7119385973 0.0000000000
+ V -0.6358200627 0.0000000000
+ V -0.0703936743 0.0000000000
+ N 0.7438340890 0.0000000000
+ N -0.6643053754 0.0000000000
+ N -0.0735473745 0.0000000000
+ W 1.8175356507 0.0000000000
+ W -1.6232096923 0.0000000000
+ W -0.1797107407 0.0000000000
+
+ E 146.5031967184 0.0000000000 11.16
+ V 0.2265837636 0.0000000000
+ V 0.3474697082 0.0000000000
+ V -0.8468727535 0.0000000000
+ N 0.2402770238 0.0000000000
+ N 0.3684685345 0.0000000000
+ N -0.8980522763 0.0000000000
+ W 2.9082771721 0.0000000000
+ W 4.4598880523 0.0000000000
+ W -10.8698904856 0.0000000000
+
+</PRE>
+
+<H2>PROGRAM NOTES</H2>
+
+The relative importance of the eigen value (E) is the ratio (percentage)
+of the eigen value to the sum of the eigen values. Note that the output
+is not sorted by relative importance.
+
+
+<P>
+
+In general, the solution to the eigen system results in complex
+numbers (with both real and imaginary parts). However, in the example
+above, since the input matrix is symmetric (i.e., inverting the rows and columns
+gives the same matrix) the eigen system has only real values (i.e., the
+imaginary part is zero).
+This fact makes it possible to use eigen vectors to perform principle component
+transformation of data sets. The covariance or correlation
+matrix of any data set is symmetric
+and thus has only real eigen values and vectors.
+
+<H2>PRINCIPLE COMPONENTS</H2>
+
+To perform principle component transformation on GRASS data layers,
+one would use
+<EM>r.covar</EM>
+to get the covariance (or correlation) matrix for a set of data layers,
+<EM>m.eigensystem</EM>
+to extract the related eigen vectors, and
+<EM>r.mapcalc</EM>
+to form the desired components.
+For example, to get the eigen vectors for 3 layers:
+
+<PRE>
+<B>(echo 3; r.covar map.1,map.2,map.3) | m.eigensystem</B>
+</PRE>
+
+Note that since m.covar only outputs the matrix, we must manually prepend the matrix
+size (3) using the echo command.
+
+
+<P>
+
+Then, using the W vector, new maps are created:
+
+
+<PRE>
+<B>r.mapcalc</B> 'pc.1 = 21.2395*map.1 + 22.1470*map.2 + 14.7696*map.3'
+<B>r.mapcalc</B> 'pc.2 = 2.9083*map.1 + 4.4599*map.2 - 10.8699*map.3'
+<B>r.mapcalc</B> 'pc.3 = 1.8175*map.1 - 1.6232*map.2 \- 0.1797*map.3'
+</PRE>
+
+<H2>NOTES</H2>
+
+The source code for this program requires a Fortran compiler.
+
+<H2>SEE ALSO</H2>
+
+<EM><A HREF="r.covar.html">r.covar</A></EM><br>
+<EM><A HREF="r.mapcalc.html">r.mapcalc</A></EM><br>
+<EM><A HREF="r.rescale.html">r.rescale</A></EM>
+
+<H2>AUTHOR</H2>
+
+This code uses routines from the EISPACK system. The interface was coded by
+Michael Shapiro, U.S.Army Construction Engineering
+Research Laboratory
+<p><i>Last changed: $Date: 2002/01/25 05:45:33 $</i>
+</body>
+</html>
Added: grass-addons/m.eigensystem/main.f
===================================================================
--- grass-addons/m.eigensystem/main.f (rev 0)
+++ grass-addons/m.eigensystem/main.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,124 @@
+ PROGRAM EIGEN
+C
+C Computes eigenvalues and eigenvectors for an NxN matrix
+C
+C AUTHOR: Michael Shapiro, U.S.Army Construction Engineering Research Laboratory
+C
+C COPYRIGHT: (C) 1999 by the GRASS Development Team
+C
+C This program is free software under the GNU General Public
+C License (>=v2). Read the file COPYING that comes with GRASS
+C for details.
+C
+C The input is from stdin and has the following format
+C
+C First line contains one positive integer N, the diminension of the matrix
+C The the matrix of real values follows N reals per line each separated by
+C white space:
+C
+C 3
+C 1.2 2.3 3.4
+C 0.1 12.3 34.5676
+C 2 12.45 2
+C
+C
+C Currently, the max value for N is 30
+C
+C The output is N sets of values. One E line and N V lines
+C
+C E real imaginary percent-importance
+C V real imaginary
+C N real imaginary
+C ...
+C
+C where E is the eigen value (and it relative importance)
+C and V are the eigenvector for this eigenvalue.
+C N are the normalized eigenvector for this eigenvalue.
+C
+ DOUBLE PRECISION A(30,30),W1,W2,WR(30),WI(30),Z(30,30),FV1(30)
+ DOUBLE PRECISION SUM1,SUM2
+ INTEGER IV1(30),N
+
+C
+C read the matrix size
+ READ(5,*,END=1000) N
+ IF(N.LT.1) GOTO 2000
+ IF(N.GT.30) GOTO 3000
+C PRINT *,"N=",N
+
+C
+C read the matrix
+ READ(5,*,END=1000) (( A(I,J),J=1,N),I=1,N)
+
+C
+C run the real-general eigen subroutine
+ CALL RG(30,N,A,WR,WI,1,Z,IV1,FV1,IERR)
+ IF (IERR.NE.0) WRITE(6,*) "? ERROR CODE",IERR
+
+
+C
+C WR and WI contain the real and imaginary parts of the eigenvalues
+C Z contains the vectors
+ SUM1 = 0.0
+ DO 1 I=1,N
+ 1 SUM1 = SUM1 + SQRT(WR(I)*WR(I)+WI(I)*WI(I))
+ DO 2 I=1,N
+ W1 = SQRT(WR(I)*WR(I)+WI(I)*WI(I))
+ W2 = SQRT(W1)
+ WRITE(6,100) "E",WR(I),WI(I), W1/SUM1 * 100.0
+
+C
+C Normalize the eigenvectors before printing
+ SUM2 = 0.0
+ DO 6 J=1,N
+ IF(WI(I))3,4,5
+ 3 SUM2 = SUM2 + Z(J,I-1) * Z(J,I-1) + Z(J,I) * Z(J,I)
+ GOTO 6
+ 4 SUM2 = SUM2 + Z(J,I) * Z(J,I)
+ GOTO 6
+ 5 SUM2 = SUM2 + Z(J,I+1) * Z(J,I+1) + Z(J,I) * Z(J,I)
+ 6 CONTINUE
+ SUM2 = SQRT(SUM2)
+
+
+ DO 16 J=1,N
+ IF(WI(I))13,14,15
+ 13 WRITE(6,100) "V",Z(J,I-1),-Z(J,I)
+ GOTO 16
+ 14 WRITE(6,100) "V",Z(J,I),0.0
+ GOTO 16
+ 15 WRITE(6,100) "V",Z(J,I),Z(J,I+1)
+ 16 CONTINUE
+
+ DO 26 J=1,N
+ IF(WI(I))23,24,25
+ 23 WRITE(6,100) "N",Z(J,I-1)/SUM2,-Z(J,I)/SUM2
+ GOTO 26
+ 24 WRITE(6,100) "N",Z(J,I)/SUM2,0.0
+ GOTO 26
+ 25 WRITE(6,100) "N",Z(J,I)/SUM2,Z(J,I+1)/SUM2
+ 26 CONTINUE
+
+ DO 36 J=1,N
+ IF(WI(I))33,34,35
+ 33 WRITE(6,100) "W",W2*Z(J,I-1)/SUM2,-W2*Z(J,I)/SUM2
+ GOTO 36
+ 34 WRITE(6,100) "W",W2*Z(J,I)/SUM2,0.0
+ GOTO 36
+ 35 WRITE(6,100) "W",W2*Z(J,I)/SUM2,W2*Z(J,I+1)/SUM2
+ 36 CONTINUE
+
+ 2 CONTINUE
+
+
+C
+ CALL EXIT(0)
+ 1000 PRINT *,"Incomplete input file"
+ CALL EXIT(1)
+ 2000 PRINT *,"N must be positive"
+ CALL EXIT(1)
+ 3000 PRINT *,"Maximum array size is 30"
+ CALL EXIT(1)
+
+ 100 FORMAT(A1,1x,F20.10,1x,F20.10,3x,F6.2)
+ END
Added: grass-addons/m.eigensystem/matrix.7x7
===================================================================
--- grass-addons/m.eigensystem/matrix.7x7 (rev 0)
+++ grass-addons/m.eigensystem/matrix.7x7 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,8 @@
+7
+1 2 3 4 5 6 7
+3 4 2 5 1 6 7
+2 1 4 2 7 6 8
+4 0 0 2 8 9 1
+5 1 2 1 2 1 0
+6 9 8 7 6 5 4
+7 0 9 8 0 5 7
Added: grass-addons/m.eigensystem/rg.f
===================================================================
--- grass-addons/m.eigensystem/rg.f (rev 0)
+++ grass-addons/m.eigensystem/rg.f 2008-01-24 08:50:35 UTC (rev 29816)
@@ -0,0 +1,70 @@
+ SUBROUTINE RG(NM,N,A,WR,WI,MATZ,Z,IV1,FV1,IERR) EIS8
+C EIS8
+ INTEGER N,NM,IS1,IS2,IERR,MATZ EIS8
+ DOUBLE PRECISION A(NM,N),WR(N),WI(N),Z(NM,N),FV1(N) EIS8
+ INTEGER IV1(N) EIS8
+C EIS8
+C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF EIS8
+C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) EIS8
+C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) EIS8
+C OF A REAL GENERAL MATRIX. EIS8
+C EIS8
+C ON INPUT EIS8
+C EIS8
+C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL EIS8
+C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM EIS8
+C DIMENSION STATEMENT. EIS8
+C EIS8
+C N IS THE ORDER OF THE MATRIX A. EIS8
+C EIS8
+C A CONTAINS THE REAL GENERAL MATRIX. EIS8
+C EIS8
+C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF EIS8
+C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO EIS8
+C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. EIS8
+C EIS8
+C ON OUTPUT EIS8
+C EIS8
+C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, EIS8
+C RESPECTIVELY, OF THE EIGENVALUES. COMPLEX CONJUGATE EIS8
+C PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY WITH THE EIS8
+C EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST. EIS8
+C EIS8
+C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS EIS8
+C IF MATZ IS NOT ZERO. IF THE J-TH EIGENVALUE IS REAL, THE EIS8
+C J-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. IF THE J-TH EIS8
+C EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY PART, THE EIS8
+C J-TH AND (J+1)-TH COLUMNS OF Z CONTAIN THE REAL AND EIS8
+C IMAGINARY PARTS OF ITS EIGENVECTOR. THE CONJUGATE OF THIS EIS8
+C VECTOR IS THE EIGENVECTOR FOR THE CONJUGATE EIGENVALUE. EIS8
+C EIS8
+C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR EIS8
+C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR HQR EIS8
+C AND HQR2. THE NORMAL COMPLETION CODE IS ZERO. EIS8
+C EIS8
+C IV1 AND FV1 ARE TEMPORARY STORAGE ARRAYS. EIS8
+C EIS8
+C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, EIS8
+C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY EIS8
+C EIS8
+C THIS VERSION DATED AUGUST 1983. EIS8
+C EIS8
+C ------------------------------------------------------------------EIS8
+C EIS8
+ IF (N .LE. NM) GO TO 10 EIS8
+ IERR = 10 * N EIS8
+ GO TO 50 EIS8
+C EIS8
+ 10 CALL BALANC(NM,N,A,IS1,IS2,FV1) EIS8
+ CALL ELMHES(NM,N,IS1,IS2,A,IV1) EIS8
+ IF (MATZ .NE. 0) GO TO 20 EIS8
+C .......... FIND EIGENVALUES ONLY .......... EIS8
+ CALL HQR(NM,N,IS1,IS2,A,WR,WI,IERR) EIS8
+ GO TO 50 EIS8
+C .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... EIS8
+ 20 CALL ELTRAN(NM,N,IS1,IS2,A,IV1,Z) EIS8
+ CALL HQR2(NM,N,IS1,IS2,A,WR,WI,Z,IERR) EIS8
+ IF (IERR .NE. 0) GO TO 50 EIS8
+ CALL BALBAK(NM,N,IS1,IS2,FV1,N,Z) EIS8
+ 50 RETURN EIS8
+ END EIS8