Mercurial > hg > octave-nkf
annotate liboctave/cruft/daspk/dorth.f @ 20830:b65888ec820e draft default tip gccjit
dmalcom gcc jit import
| author | Stefan Mahr <dac922@gmx.de> |
|---|---|
| date | Fri, 27 Feb 2015 16:59:36 +0100 |
| parents | 446c46af4b42 |
| children |
| rev | line source |
|---|---|
| 3911 | 1 C Work performed under the auspices of the U.S. Department of Energy |
|
19790
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
15271
diff
changeset
|
2 C by Lawrence Livermore National Laboratory under contract number |
| 3911 | 3 C W-7405-Eng-48. |
| 4 C | |
| 5 SUBROUTINE DORTH (VNEW, V, HES, N, LL, LDHES, KMP, SNORMW) | |
| 6 C | |
| 7 C***BEGIN PROLOGUE DORTH | |
| 8 C***DATE WRITTEN 890101 (YYMMDD) | |
| 9 C***REVISION DATE 900926 (YYMMDD) | |
| 10 C | |
| 11 C | |
| 12 C----------------------------------------------------------------------- | |
| 13 C***DESCRIPTION | |
| 14 C | |
| 15 C This routine orthogonalizes the vector VNEW against the previous | |
| 16 C KMP vectors in the V array. It uses a modified Gram-Schmidt | |
| 17 C orthogonalization procedure with conditional reorthogonalization. | |
| 18 C | |
| 19 C On entry | |
| 20 C | |
| 21 C VNEW = The vector of length N containing a scaled product | |
| 22 C OF The Jacobian and the vector V(*,LL). | |
| 23 C | |
| 24 C V = The N x LL array containing the previous LL | |
| 25 C orthogonal vectors V(*,1) to V(*,LL). | |
| 26 C | |
| 27 C HES = An LL x LL upper Hessenberg matrix containing, | |
| 28 C in HES(I,K), K.LT.LL, scaled inner products of | |
| 29 C A*V(*,K) and V(*,I). | |
| 30 C | |
| 31 C LDHES = The leading dimension of the HES array. | |
| 32 C | |
| 33 C N = The order of the matrix A, and the length of VNEW. | |
| 34 C | |
| 35 C LL = The current order of the matrix HES. | |
| 36 C | |
| 37 C KMP = The number of previous vectors the new vector VNEW | |
| 38 C must be made orthogonal to (KMP .LE. MAXL). | |
| 39 C | |
| 40 C | |
| 41 C On return | |
| 42 C | |
| 43 C VNEW = The new vector orthogonal to V(*,I0), | |
| 44 C where I0 = MAX(1, LL-KMP+1). | |
| 45 C | |
| 46 C HES = Upper Hessenberg matrix with column LL filled in with | |
| 47 C scaled inner products of A*V(*,LL) and V(*,I). | |
| 48 C | |
| 49 C SNORMW = L-2 norm of VNEW. | |
| 50 C | |
| 51 C----------------------------------------------------------------------- | |
| 52 C***ROUTINES CALLED | |
|
19790
446c46af4b42
strip trailing whitespace from most source files
John W. Eaton <jwe@octave.org>
parents:
15271
diff
changeset
|
53 C DDOT, DNRM2, DAXPY |
| 3911 | 54 C |
| 55 C***END PROLOGUE DORTH | |
| 56 C | |
| 57 INTEGER N, LL, LDHES, KMP | |
| 58 DOUBLE PRECISION VNEW, V, HES, SNORMW | |
| 59 DIMENSION VNEW(*), V(N,*), HES(LDHES,*) | |
| 60 INTEGER I, I0 | |
| 61 DOUBLE PRECISION ARG, DDOT, DNRM2, SUMDSQ, TEM, VNRM | |
| 62 C | |
| 63 C----------------------------------------------------------------------- | |
| 64 C Get norm of unaltered VNEW for later use. | |
| 65 C----------------------------------------------------------------------- | |
| 66 VNRM = DNRM2 (N, VNEW, 1) | |
| 67 C----------------------------------------------------------------------- | |
| 68 C Do Modified Gram-Schmidt on VNEW = A*V(LL). | |
| 69 C Scaled inner products give new column of HES. | |
| 70 C Projections of earlier vectors are subtracted from VNEW. | |
| 71 C----------------------------------------------------------------------- | |
| 72 I0 = MAX0(1,LL-KMP+1) | |
| 73 DO 10 I = I0,LL | |
| 74 HES(I,LL) = DDOT (N, V(1,I), 1, VNEW, 1) | |
| 75 TEM = -HES(I,LL) | |
| 76 CALL DAXPY (N, TEM, V(1,I), 1, VNEW, 1) | |
| 77 10 CONTINUE | |
| 78 C----------------------------------------------------------------------- | |
| 79 C Compute SNORMW = norm of VNEW. | |
| 80 C If VNEW is small compared to its input value (in norm), then | |
| 81 C Reorthogonalize VNEW to V(*,1) through V(*,LL). | |
| 82 C Correct if relative correction exceeds 1000*(unit roundoff). | |
| 83 C Finally, correct SNORMW using the dot products involved. | |
| 84 C----------------------------------------------------------------------- | |
| 85 SNORMW = DNRM2 (N, VNEW, 1) | |
| 86 IF (VNRM + 0.001D0*SNORMW .NE. VNRM) RETURN | |
| 87 SUMDSQ = 0.0D0 | |
| 88 DO 30 I = I0,LL | |
| 89 TEM = -DDOT (N, V(1,I), 1, VNEW, 1) | |
| 90 IF (HES(I,LL) + 0.001D0*TEM .EQ. HES(I,LL)) GO TO 30 | |
| 91 HES(I,LL) = HES(I,LL) - TEM | |
| 92 CALL DAXPY (N, TEM, V(1,I), 1, VNEW, 1) | |
| 93 SUMDSQ = SUMDSQ + TEM**2 | |
| 94 30 CONTINUE | |
| 95 IF (SUMDSQ .EQ. 0.0D0) RETURN | |
| 96 ARG = MAX(0.0D0,SNORMW**2 - SUMDSQ) | |
| 97 SNORMW = SQRT(ARG) | |
| 98 RETURN | |
| 99 C | |
| 100 C------END OF SUBROUTINE DORTH------------------------------------------ | |
| 101 END |