-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathderi21.f
executable file
·63 lines (63 loc) · 2.02 KB
/
deri21.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
SUBROUTINE DERI21 (A,NVAR,MINEAR,NFIRST,VNERT,PNERT
1 ,B,NCUT)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(MINEAR,NVAR),VNERT(NVAR),PNERT(NVAR),B(MINEAR,*)
************************************************************************
*
* LEAST-SQUARE ANALYSIS OF A SET OF NVAR POINTS {A} :
*
* PRODUCE A SUBSET OF NCUT ORTHONORMALIZED VECTORS B, OPTIMUM IN A
* LEAST-SQUARE SENSE WITH RESPECT TO THE INITIAL SPACE {A}.
* EACH NEW HIERARCHIZED VECTOR B EXTRACTS A MAXIMUM PERCENTAGE FROM
* THE REMAINING DISPERSION OF THE SET {A} OUT OF THE PREVIOUS
* {B} SUBSPACE.
* INPUT
* A(MINEAR,NVAR): ORIGINAL SET {A}.
* NFIRST : MAXIMUM ALLOWED SIZE OF THE BASIS B.
* OUTPUT
* VNERT(NVAR) : LOWEST EIGENVECTOR OF A'* A.
* PNERT(NVAR) : SQUARE ROOT OF THE ASSOCIATED EIGENVALUES
* IN DECREASING ORDER.
* B(MINEAR,NCUT): OPTIMUM ORTHONORMALIZED SUBSET {B}.
*
************************************************************************
DIMENSION WORK(4)
C
C VNERT = A' * A
CUTOFF=0.85D0
SUM2=0.D0
CALL MTXMC(A,NVAR,A,MINEAR,WORK)
DO 10 I=1,(NVAR*(NVAR+1))/2
10 WORK(I)=-WORK(I)
C DIAGONALIZE IN DECREASING ORDER OF EIGENVALUES
IF(ABS(WORK(1)).LT.1.D-28 .AND. NVAR.EQ.1)THEN
PNERT(1)=SQRT(-WORK(1))
WORK(1)=1.D15
VNERT(1)=1.D0
NCUT=1
GOTO 50
ELSE
CALL HQRII(WORK,NVAR,NVAR,PNERT, VNERT)
C FIND NCUT ACCORDING TO CUTOFF, BUILD WORK = VNERT * (PNERT)**-0.5
SUM=0.D0
DO 20 I=1,NVAR
20 SUM=SUM-PNERT(I)
L=1
DO 40 I=1,NFIRST
SUM2=SUM2-PNERT(I)/SUM
PNERT(I)=SQRT(-PNERT(I))
DO 30 J=1,NVAR
WORK(L)=VNERT(L)/PNERT(I)
30 L=L+1
IF(SUM2.GE.CUTOFF) THEN
NCUT=I
GO TO 50
ENDIF
40 CONTINUE
NCUT=NFIRST
C ORTHONORMALIZED BASIS
C B(MINEAR,NCUT) = A(MINEAR,NVAR)*WORK(NVAR,NCUT)
ENDIF
50 CALL MXM (A,MINEAR,WORK,NVAR,B,NCUT)
RETURN
END