#include /* nagroutines.f -- translated by f2c (version 19970805). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Table of constant values */ static doublereal c_b15 = 1.; static integer c__1 = 1; static doublereal c_b24 = -1.; static integer c__0 = 0; static integer c_n1 = -1; static doublereal c_b114 = 0.; static logical c_false = FALSE_; static logical c_true = TRUE_; doublereal a02abf_(xxr, xxi) doublereal *xxr, *xxi; { /* Initialized data */ static doublereal zero = 0.; static doublereal one = 1.; /* System generated locals */ doublereal ret_val, d__1; /* Builtin functions */ double sqrt(); /* Local variables */ static doublereal h__, xi, xr; /* NAG COPYRIGHT 1975 */ /* MARK 4.5 REVISED */ /* MARK 5C REVISED */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* RETURNS THE ABSOLUTE VALUE OF A COMPLEX NUMBER VIA ROUTINE */ /* NAME */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Intrinsic Functions .. */ /* .. Data statements .. */ /* .. Executable Statements .. */ xr = abs(*xxr); xi = abs(*xxi); if (xi <= xr) { goto L20; } h__ = xr; xr = xi; xi = h__; L20: if (xi != zero) { goto L40; } ret_val = xr; return ret_val; L40: /* Computing 2nd power */ d__1 = xi / xr; h__ = xr * sqrt(one + d__1 * d__1); ret_val = h__; return ret_val; } /* a02abf_ */ /* Subroutine */ int a02acf_(xxr, xxi, yyr, yyi, zr, zi) doublereal *xxr, *xxi, *yyr, *yyi, *zr, *zi; { /* Initialized data */ static doublereal one = 1.; static doublereal a, h__; /* MARK 2A RELEASE. NAG COPYRIGHT 1973 */ /* MARK 4.5 REVISED */ /* MARK 5C REVISED */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* DIVIDES ONE COMPLEX NUMBER BY A SECOND */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Intrinsic Functions .. */ /* .. Data statements .. */ /* .. Executable Statements .. */ if (abs(*yyr) <= abs(*yyi)) { goto L20; } h__ = *yyi / *yyr; a = one / (h__ * *yyi + *yyr); *zr = (*xxr + h__ * *xxi) * a; *zi = (*xxi - h__ * *xxr) * a; return 0; L20: h__ = *yyr / *yyi; a = one / (h__ * *yyr + *yyi); *zr = (h__ * *xxr + *xxi) * a; *zi = (h__ * *xxi - *xxr) * a; return 0; } /* a02acf_ */ /* Subroutine */ int f01akf_(n, k, l, a, ia, intger) integer *n, *k, *l; doublereal *a; integer *ia, *intger; { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; doublereal d__1; /* Local variables */ static integer i__, j, m; static doublereal x, y; extern /* Subroutine */ int dgemv_(); static integer k1; extern /* Subroutine */ int dtrsv_(); /* MARK 2 RELEASE. NAG COPYRIGHT 1972 */ /* MARK 4 REVISED. */ /* MARK 4.5 REVISED */ /* MARK 8 REVISED. IER-248 (JUN 1980). */ /* MARK 11 REVISED. VECTORISATION (JAN 1984). */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 12 REVISED. EXTENDED BLAS (JUNE 1986) */ /* DIRHES */ /* AUGUST 1ST, 1971 . */ /* GIVEN THE UNSYMMETRIC MATRIX, A, STORED IN THE ARRAY A(N,N), */ /* THIS SUBROUTINE REDUCES THE SUB-MATRIX OF ORDER L - K + 1, */ /* WHICH STARTS AT THE ELEMENT A(K,K) AND FINISHES AT THE */ /* ELEMENT A(L,L), TO HESSENBERG FORM, H, BY THE DIRECT */ /* METHOD(AN = NH). THE MATRIX H IS OVERWRITTEN ON A WITH */ /* DETAILS OF THE TRANSFORMATIONS (N) STORED IN THE REMAINING */ /* TRIANGLE UNDER H AND IN ELEMENTS K TO L OF THE ARRAY */ /* INTGER(N). */ /* 1ST AUGUST 1971 */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --intger; a_dim1 = *ia; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ k1 = *k + 1; if (k1 > *l) { return 0; } i__1 = *n; for (j = k1; j <= i__1; ++j) { m = j; x = 0.; if (j > *l) { goto L120; } i__2 = *l; for (i__ = j; i__ <= i__2; ++i__) { if ((d__1 = a[i__ + (j - 1) * a_dim1], abs(d__1)) <= abs(x)) { goto L20; } x = a[i__ + (j - 1) * a_dim1]; m = i__; L20: ; } intger[j] = m; if (m == j) { goto L80; } /* INTERCHANGE ROWS AND COLUMNS OF A. */ i__2 = *n; for (i__ = *k; i__ <= i__2; ++i__) { y = a[m + i__ * a_dim1]; a[m + i__ * a_dim1] = a[j + i__ * a_dim1]; a[j + i__ * a_dim1] = y; /* L40: */ } i__2 = *l; for (i__ = 1; i__ <= i__2; ++i__) { y = a[i__ + m * a_dim1]; a[i__ + m * a_dim1] = a[i__ + j * a_dim1]; a[i__ + j * a_dim1] = y; /* L60: */ } L80: if (x != 0. && j < *l) { i__2 = *l; for (i__ = j + 1; i__ <= i__2; ++i__) { a[i__ + (j - 1) * a_dim1] /= x; /* L100: */ } i__2 = *l - j; dgemv_("N", l, &i__2, &c_b15, &a[(j + 1) * a_dim1 + 1], ia, &a[j + 1 + (j - 1) * a_dim1], &c__1, &c_b15, &a[j * a_dim1 + 1] , &c__1, 1L); } L120: i__2 = j - *k; dtrsv_("L", "N", "U", &i__2, &a[*k + 1 + *k * a_dim1], ia, &a[*k + 1 + j * a_dim1], &c__1, 1L, 1L, 1L); if (j < *l) { i__2 = *l - j; i__3 = j - *k; dgemv_("N", &i__2, &i__3, &c_b24, &a[j + 1 + *k * a_dim1], ia, &a[ *k + 1 + j * a_dim1], &c__1, &c_b15, &a[j + 1 + j * a_dim1], &c__1, 1L); } /* L140: */ } return 0; } /* f01akf_ */ /* Subroutine */ int f01apf_(n, low, iupp, intger, h__, ih, v, iv) integer *n, *low, *iupp, *intger; doublereal *h__; integer *ih; doublereal *v; integer *iv; { /* System generated locals */ integer h_dim1, h_offset, v_dim1, v_offset, i__1, i__2; /* Local variables */ static integer i__, j, m; static doublereal x; static integer i1, ii, low1; /* MARK 2 RELEASE. NAG COPYRIGHT 1972 */ /* MARK 4 REVISED. */ /* MARK 4.5 REVISED */ /* MARK 5C REVISED */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* DIRTRANS */ /* FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS IN THE ARRAY */ /* V(N,N) FROM THE INFORMATION LEFT BY SUBROUTINE F01AKF */ /* BELOW THE UPPER HESSENBERG MATRIX, H, IN THE ARRAY H(N,N) */ /* AND IN THE INTEGER ARRAY INTGER(N). */ /* 1ST AUGUST 1971 */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --intger; h_dim1 = *ih; h_offset = h_dim1 + 1; h__ -= h_offset; v_dim1 = *iv; v_offset = v_dim1 + 1; v -= v_offset; /* Function Body */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *n; for (j = 1; j <= i__2; ++j) { v[i__ + j * v_dim1] = 0.; /* L20: */ } v[i__ + i__ * v_dim1] = 1.; /* L40: */ } low1 = *low + 1; if (low1 > *iupp) { return 0; } i__1 = *iupp; for (ii = low1; ii <= i__1; ++ii) { i__ = low1 + *iupp - ii; i1 = i__ - 1; if (low1 > i1) { goto L80; } i__2 = i1; for (j = low1; j <= i__2; ++j) { v[i__ + j * v_dim1] = h__[i__ + (j - 1) * h_dim1]; /* L60: */ } L80: m = intger[i__]; if (m == i__) { goto L120; } i__2 = *iupp; for (j = low1; j <= i__2; ++j) { x = v[m + j * v_dim1]; v[m + j * v_dim1] = v[i__ + j * v_dim1]; v[i__ + j * v_dim1] = x; /* L100: */ } L120: ; } return 0; } /* f01apf_ */ /* Subroutine */ int f01atf_(n, ib, a, ia, low, lhi, d__) integer *n, *ib; doublereal *a; integer *ia, *low, *lhi; doublereal *d__; { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2; doublereal d__1; /* Local variables */ static doublereal c__, f, g; static integer i__, j, k, l; static doublereal r__, s; extern /* Subroutine */ int f01atz_(); static doublereal b2; static integer jj; static logical noconv; /* MARK 2 RELEASE. NAG COPYRIGHT 1972 */ /* MARK 4 REVISED. */ /* MARK 4.5 REVISED */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* BALANCE */ /* REDUCE THE NORM OF A(N,N) BY EXACT DIAGONAL SIMILARITY */ /* TRANSFORMATIONS STORED IN D(N). */ /* DECEMBER 1ST.,1971 */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --d__; a_dim1 = *ia; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ b2 = (doublereal) (*ib * *ib); l = 1; k = *n; L20: if (k < 1) { goto L100; } /* SEARCH FOR ROWS ISOLATING AN EIGENVALUE AND PUSH THEM DOWN */ j = k + 1; i__1 = k; for (jj = 1; jj <= i__1; ++jj) { --j; r__ = 0.; i__2 = k; for (i__ = 1; i__ <= i__2; ++i__) { if (i__ == j) { goto L40; } r__ += (d__1 = a[j + i__ * a_dim1], abs(d__1)); L40: ; } if (r__ == 0.) { goto L80; } /* L60: */ } goto L100; L80: f01atz_(&k, &a[a_offset], ia, &d__[1], &k, &l, n, &j); --k; goto L20; /* SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE AND PUSH THEM */ /* LEFT. */ L100: if (l > k) { goto L180; } i__1 = k; for (j = l; j <= i__1; ++j) { c__ = 0.; i__2 = k; for (i__ = l; i__ <= i__2; ++i__) { if (i__ == j) { goto L120; } c__ += (d__1 = a[i__ + j * a_dim1], abs(d__1)); L120: ; } if (c__ == 0.) { goto L160; } /* L140: */ } goto L180; L160: f01atz_(&l, &a[a_offset], ia, &d__[1], &k, &l, n, &j); ++l; goto L100; /* NOW BALANCE THE SUBMATRIX IN ROWS L THROUGH K. */ L180: *low = l; *lhi = k; if (l > k) { goto L220; } i__1 = k; for (i__ = l; i__ <= i__1; ++i__) { d__[i__] = 1.; /* L200: */ } L220: noconv = FALSE_; if (l > k) { goto L420; } i__1 = k; for (i__ = l; i__ <= i__1; ++i__) { c__ = 0.; r__ = 0.; i__2 = k; for (j = l; j <= i__2; ++j) { if (j == i__) { goto L240; } c__ += (d__1 = a[j + i__ * a_dim1], abs(d__1)); r__ += (d__1 = a[i__ + j * a_dim1], abs(d__1)); L240: ; } g = r__ / (doublereal) (*ib); f = 1.; s = c__ + r__; L260: if (c__ >= g) { goto L280; } f *= (doublereal) (*ib); c__ *= b2; goto L260; L280: g = r__ * (doublereal) (*ib); L300: if (c__ < g) { goto L320; } f /= (doublereal) (*ib); c__ /= b2; goto L300; L320: if ((c__ + r__) / f >= s * .95) { goto L400; } g = 1. / f; d__[i__] *= f; noconv = TRUE_; if (l > *n) { goto L360; } i__2 = *n; for (j = l; j <= i__2; ++j) { a[i__ + j * a_dim1] *= g; /* L340: */ } L360: i__2 = k; for (j = 1; j <= i__2; ++j) { a[j + i__ * a_dim1] *= f; /* L380: */ } L400: ; } L420: if (noconv) { goto L220; } return 0; } /* f01atf_ */ /* Subroutine */ int f01atz_(m, a, ia, d__, k, l, n, j) integer *m; doublereal *a; integer *ia; doublereal *d__; integer *k, *l, *n, *j; { /* System generated locals */ integer a_dim1, a_offset, i__1; /* Local variables */ static doublereal f; static integer i__; /* MARK 2 RELEASE. NAG COPYRIGHT 1972 */ /* MARK 4 REVISED. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* NAG COPYRIGHT 1975 */ /* MARK 4.5 REVISED */ /* AUXILIARY ROUTINE CALLED BY F01ATF. */ /* INTERCHANGES ELEMENTS 1 TO K OF COLUMNS J AND M, */ /* AND ELEMENTS L TO N OF ROWS J AND M. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --d__; a_dim1 = *ia; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ d__[*m] = (doublereal) (*j); if (*j == *m) { goto L60; } i__1 = *k; for (i__ = 1; i__ <= i__1; ++i__) { f = a[i__ + *j * a_dim1]; a[i__ + *j * a_dim1] = a[i__ + *m * a_dim1]; a[i__ + *m * a_dim1] = f; /* L20: */ } if (*l > *n) { goto L60; } i__1 = *n; for (i__ = *l; i__ <= i__1; ++i__) { f = a[*j + i__ * a_dim1]; a[*j + i__ * a_dim1] = a[*m + i__ * a_dim1]; a[*m + i__ * a_dim1] = f; /* L40: */ } L60: return 0; } /* f01atz_ */ /* Subroutine */ int f01auf_(n, low, lhi, m, d__, z__, iz) integer *n, *low, *lhi, *m; doublereal *d__, *z__; integer *iz; { /* System generated locals */ integer z_dim1, z_offset, i__1, i__2; /* Local variables */ static integer i__, j, k; static doublereal s; static integer ii, lhi1, low1; /* MARK 2 RELEASE. NAG COPYRIGHT 1972 */ /* MARK 4 REVISED. */ /* MARK 4.5 REVISED */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* BALBAK */ /* BACKWARD TRANSFORMATION OF A SET OF RIGHT-HAND EIGENVECTORS */ /* OF A BALANCED MATRIX INTO THE EIGENVECTORS OF THE ORIGINAL */ /* MATRIX FROM WHICH THE BALANCED MATRIX WAS DERIVED BY A CALL */ /* OF SUBROUTINE F01ATF. */ /* DECEMBER 1ST.,1971 */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --d__; z_dim1 = *iz; z_offset = z_dim1 + 1; z__ -= z_offset; /* Function Body */ if (*low > *lhi) { goto L60; } i__1 = *lhi; for (i__ = *low; i__ <= i__1; ++i__) { s = d__[i__]; /* LEFT-HAND EIGENVECTORS ARE BACK TRANSFORMED IF THE */ /* FOREGOING STATEMENT IS REPLACED BY S=1/D(I) */ i__2 = *m; for (j = 1; j <= i__2; ++j) { z__[i__ + j * z_dim1] *= s; /* L20: */ } /* L40: */ } L60: i__ = *low; low1 = *low - 1; if (low1 < 1) { goto L120; } i__1 = low1; for (ii = 1; ii <= i__1; ++ii) { --i__; k = (integer) d__[i__]; if (k == i__) { goto L100; } i__2 = *m; for (j = 1; j <= i__2; ++j) { s = z__[i__ + j * z_dim1]; z__[i__ + j * z_dim1] = z__[k + j * z_dim1]; z__[k + j * z_dim1] = s; /* L80: */ } L100: ; } L120: lhi1 = *lhi + 1; if (lhi1 > *n) { return 0; } i__1 = *n; for (i__ = lhi1; i__ <= i__1; ++i__) { k = (integer) d__[i__]; if (k == i__) { goto L160; } i__2 = *m; for (j = 1; j <= i__2; ++j) { s = z__[i__ + j * z_dim1]; z__[i__ + j * z_dim1] = z__[k + j * z_dim1]; z__[k + j * z_dim1] = s; /* L140: */ } L160: ; } return 0; } /* f01auf_ */ /* Subroutine */ int f01lzf_(n, a, nra, c__, nrc, wantb, b, wantq, wanty, y, nry, ly, wantz, z__, nrz, ncz, d__, e, work1, work2, ifail) integer *n; doublereal *a; integer *nra; doublereal *c__; integer *nrc; logical *wantb; doublereal *b; logical *wantq, *wanty; doublereal *y; integer *nry, *ly; logical *wantz; doublereal *z__; integer *nrz, *ncz; doublereal *d__, *e, *work1, *work2; integer *ifail; { /* System generated locals */ integer a_dim1, a_offset, c_dim1, c_offset, y_dim1, y_offset, z_dim1, z_offset, i__1, i__2, i__3; /* Builtin functions */ double sqrt(); /* Local variables */ static integer ierr, i__, j, k; extern integer p01abf_(); static doublereal t, w, x; static char p01rec[1*1]; extern doublereal x02ajf_(), x02amf_(); static doublereal small; extern /* Subroutine */ int f01lzw_(), f01lzx_(), f01lzy_(); extern doublereal f01lzz_(); static integer jj; static doublereal cs, sn; static integer jp1, kp1, kp2, nm2; static doublereal sqteps, rsqtps, big, eps; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988). */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (BIDIAG) */ /* F01LZF RETURNS ALL OR PART OF THE FACTORIZATION OF THE */ /* N*N UPPER TRIANGULAR MATRIX A GIVEN BY */ /* A = Q*C*(P**T) , */ /* WHERE Q AND P ARE N*N ORTHOGONAL MATRICES AND C IS AN */ /* N*N UPPER BIDIAGONAL MATRIX. */ /* IF WANTB IS .TRUE. THEN B RETURNS (Q**T)*B. */ /* IF WANTY IS .TRUE. THEN Y RETURNS Y*Q. */ /* IF WANTZ IS .TRUE. THEN Z RETURNS (P**T)*Z. */ /* INPUT PARAMETERS. */ /* N - ORDER OF THE MATRIX A. */ /* A - THE N*N UPPER TRIANGULAR MATRIX TO BE FACTORIZED. THE */ /* STRICTLY LOWER TRIANGULAR PART OF A IS NOT REFERENCED. */ /* NRA - ROW DIMENSION OF A AS DECLARED IN THE CALLING PROGRAM. */ /* NRA MUST BE AT LEAST N. */ /* NRC - ROW DIMENSION OF C AS DECLARED IN THE CALLING PROGRAM. */ /* NRC MUST BE AT LEAST N. */ /* WANTB - MUST BE .TRUE. IF (Q**T)*B IS REQUIRED. */ /* IF WANTB IS .FALSE. THEN B IS NOT REFERENCED. */ /* B - AN N ELEMENT REAL VECTOR. */ /* WANTQ - MUST BE .TRUE. IF DETAILS OF Q ARE TO BE */ /* STORED BELOW THE BIDIAGONAL PART OF C. */ /* IF WANTQ IS .FALSE. THEN THE LOWER TRIANGULAR */ /* PART OF C IS NOT REFERENCED. */ /* WANTY - MUST BE .TRUE. IF Y*Q IS REQUIRED. */ /* IF WANTY IS .FALSE. THEN Y IS NOT REFERENCED. */ /* Y - AN LY*N REAL MATRIX. */ /* NRY - IF WANTY IS .TRUE. THEN NRY MUST BE THE ROW */ /* DIMENSION OF Y AS DECLARED IN THE CALLING */ /* PROGRAM AND MUST BE AT LEAST LY. */ /* LY - IF WANTY IS .TRUE. THEN LY MUST BE THE NUMBER */ /* OF ROWS OF Y AND MUST BE AT LEAST 1. */ /* WANTZ - MUST BE .TRUE. IF (P**T)*Z IS REQUIRED. */ /* IF WANTZ IS .FALSE. THEN Z IS NOT REFERENCED. */ /* Z - AN N*NCZ REAL MATRIX. */ /* NRZ - IF WANTZ IS .TRUE. THEN NRZ MUST BE THE ROW */ /* DIMENSION OF Z AS DECLARED IN THE CALLING */ /* PROGRAM AND MUST BE AT LEAST N. */ /* NCZ - IF WANTZ IS .TRUE. THEN NCZ MUST BE THE */ /* NUMBER OF COLUMNS OF Z AND MUST BE AT LEAST */ /* 1. */ /* IFAIL - THE USUAL FAILURE PARAMETER. IF IN DOUBT SET */ /* IFAIL TO ZERO BEFORE CALLING THIS ROUTINE. */ /* OUTPUT PARAMETERS. */ /* C - N*N MATRIX CONTAINING THE UPPER BIDIAGONAL MATRIX B. */ /* DETAILS OF P ARE STORED ABOVE THE BIDIAGONAL */ /* PART OF C. UNLESS WANTQ IS .TRUE. THE */ /* STRICTLY LOWER TRIANGULAR PART OF C IS NOT */ /* REFERENCED. */ /* THE ROUTINE MAY BE CALLED WITH C=A. */ /* B - IF WANTB IS .TRUE. THEN B WILL RETURN THE N ELEMENT */ /* VECTOR (Q**T)*B. */ /* Y - IF WANTY IS .TRUE. THEN Y WILL RETURN THE */ /* LY*N MATRIX Y*Q. */ /* Z - IF WANTZ IS .TRUE. THEN Z WILL RETURN THE N*NCZ MATRIX */ /* (P**T)*Z. */ /* D - N ELEMENT VECTOR CONTAINING THE DIAGONAL ELEMENTS OF C */ /* SUCH THAT D(I)=C(I,I), I=1,2,...,N. */ /* E - N ELEMENT VECTOR CONTAINING THE */ /* SUPER-DIAGONAL ELEMENTS OF C SUCH THAT */ /* E(I)=C(I-1,I), I=2,3,...,N. E(1) IS NOT */ /* REFERENCED. */ /* IFAIL - ON NORMAL RETURN IFAIL WILL BE ZERO. */ /* IF AN INPUT PARAMETER IS INCORRECTLY SUPPLIED */ /* THEN IFAIL IS SET TO UNITY. NO OTHER FAILURE */ /* IS POSSIBLE. */ /* WORKSPACE PARAMETERS. */ /* WORK1 */ /* WORK2 - N ELEMENT REAL VECTORS. */ /* IF WANTZ IS .FALSE. THEN WORK1 AND WORK2 ARE NOT */ /* REFERENCED. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --work2; --work1; --e; --d__; --b; a_dim1 = *nra; a_offset = a_dim1 + 1; a -= a_offset; c_dim1 = *nrc; c_offset = c_dim1 + 1; c__ -= c_offset; y_dim1 = *nry; y_offset = y_dim1 + 1; y -= y_offset; z_dim1 = *nrz; z_offset = z_dim1 + 1; z__ -= z_offset; /* Function Body */ ierr = *ifail; if (ierr == 0) { *ifail = 1; } if (*nra < *n || *nrc < *n || *n < 1) { goto L220; } if (*wanty && (*nry < *ly || *ly < 1)) { goto L220; } if (*wantz && (*nrz < *n || *ncz < 1)) { goto L220; } small = x02amf_(); big = 1. / small; eps = x02ajf_(); sqteps = sqrt(eps); rsqtps = 1. / sqteps; d__[1] = a[a_dim1 + 1]; i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = j; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] = a[i__ + j * a_dim1]; /* L20: */ } /* L40: */ } *ifail = 0; if (*n == 1) { return 0; } if (*n == 2) { goto L200; } /* START MAIN LOOP. K(TH) STEP PUTS ZEROS INTO K(TH) ROW OF C. */ nm2 = *n - 2; i__1 = nm2; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; /* SET UP PLANE ROTATION P(J,J+1) TO ANNIHILATE C(K,J+1). */ /* THIS ROTATION INTRODUCES AN UNWANTED ELEMENT IN C(J+1,J) */ /* WHICH IS STORED IN X. */ /* J GOES N-1,N-2,...,K+1. */ j = *n; i__2 = nm2; for (jj = k; jj <= i__2; ++jj) { jp1 = j; --j; w = c__[k + jp1 * c_dim1]; t = f01lzz_(&c__[k + j * c_dim1], &w, &small, &big); c__[k + jp1 * c_dim1] = t; x = 0.; f01lzw_(&t, &cs, &sn, &sqteps, &rsqtps, &big); if (! (*wantz)) { goto L60; } work1[j] = cs; work2[j] = sn; L60: if (t == 0.) { goto L80; } c__[k + j * c_dim1] = cs * c__[k + j * c_dim1] + sn * w; /* NOW APPLY THE TRANSFORMATION P(J,J+1). */ i__3 = j - k; f01lzy_(&i__3, &cs, &sn, &c__[kp1 + j * c_dim1], &c__[kp1 + jp1 * c_dim1]); x = sn * c__[jp1 + jp1 * c_dim1]; c__[jp1 + jp1 * c_dim1] = cs * c__[jp1 + jp1 * c_dim1]; /* NOW SET UP PLANE ROTATION Q(J,J+1)**T TO ANNIHILATE */ /* X=C(J+1,J). */ L80: t = f01lzz_(&c__[j + j * c_dim1], &x, &small, &big); if (*wantq) { c__[jp1 + k * c_dim1] = t; } f01lzw_(&t, &d__[j], &e[j], &sqteps, &rsqtps, &big); c__[j + j * c_dim1] = d__[j] * c__[j + j * c_dim1] + e[j] * x; if (*wanty) { f01lzy_(ly, &d__[j], &e[j], &y[j * y_dim1 + 1], &y[jp1 * y_dim1 + 1]); } /* L100: */ } /* NOW APPLY THE TRANSFORMATIONS Q(J,J+1)**T AND FORM */ /* (P(J,J+1)**T)*Z, J=N-1,N-2,...,K+1 COLUMN BY COLUMN */ kp2 = kp1 + 1; i__2 = *n; for (j = kp2; j <= i__2; ++j) { i__3 = j - k; f01lzx_(&i__3, &d__[k], &e[k], &c__[kp1 + j * c_dim1]); /* L120: */ } if (*wantb) { i__2 = *n - k; f01lzx_(&i__2, &d__[k], &e[k], &b[kp1]); } if (! (*wantz)) { goto L160; } i__2 = *ncz; for (j = 1; j <= i__2; ++j) { i__3 = *n - k; f01lzx_(&i__3, &work1[k], &work2[k], &z__[kp1 + j * z_dim1]); /* L140: */ } L160: d__[kp1] = c__[kp1 + kp1 * c_dim1]; e[kp1] = c__[k + kp1 * c_dim1]; /* L180: */ } L200: d__[*n] = c__[*n + *n * c_dim1]; e[*n] = c__[*n - 1 + *n * c_dim1]; return 0; L220: *ifail = p01abf_(&ierr, ifail, "F01LZF", &c__0, p01rec, 6L, 1L); return 0; } /* f01lzf_ */ /* Subroutine */ int f01lzw_(t, c__, s, sqteps, rsqtps, big) doublereal *t, *c__, *s, *sqteps, *rsqtps, *big; { /* Builtin functions */ double sqrt(), d_sign(); /* Local variables */ static doublereal abst, tt; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988). */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (COSSIN) */ /* F01LZW RETURNS THE VALUES */ /* C = COS(THETA) AND S = SIN(THETA) */ /* FOR A GIVEN VALUE OF */ /* T = TAN(THETA) . */ /* C IS ALWAYS NON-NEGATIVE AND S HAS THE SAME SIGN AS T. */ /* SQTEPS, RSQTPS AND BIG MUST BE SUCH THAT */ /* SQTEPS = SQRT(X02AJF) , RSQTPS = 1.0/SQTEPS AND BIG = */ /* 1.0/X02AMF , */ /* WHERE X02AJF AND X02AMF ARE THE NUMBERS RETURNED FROM */ /* ROUTINES X02AJF AND X02AMF RESPECTIVELY. */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ if (*t != 0.) { goto L20; } *c__ = 1.; *s = 0.; return 0; L20: abst = abs(*t); if (abst < *sqteps) { goto L60; } if (abst > *rsqtps) { goto L80; } tt = abst * abst; if (abst > 1.) { goto L40; } tt *= .25; *c__ = .5 / sqrt(tt + .25); *s = *c__ * *t; return 0; L40: tt = .25 / tt; *s = .5 / sqrt(tt + .25); *c__ = *s / abst; *s = d_sign(s, t); return 0; L60: *c__ = 1.; *s = *t; return 0; L80: *c__ = 0.; if (abst < *big) { *c__ = 1. / abst; } *s = d_sign(&c_b15, t); return 0; } /* f01lzw_ */ /* Subroutine */ int f01lzx_(n, c__, s, x) integer *n; doublereal *c__, *s, *x; { /* System generated locals */ integer i__1; /* Local variables */ static integer i__; static doublereal w; static integer ii, im1; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (PLROT6) */ /* F01LZX RETURNS THE N ELEMENT VECTOR */ /* Y = R(1,2)*R(2,3)*...*R(N-1,N)*X , */ /* WHERE X IS AN N ELEMENT VECTOR AND R(J-1,J) IS A PLANE */ /* ROTATION FOR THE (J-1,J)-PLANE. */ /* Y IS OVERWRITTEN ON X. */ /* THE N ELEMENT VECTORS C AND S MUST BE SUCH THAT THE */ /* NON-IDENTITY PART OF R(J-1,J) IS GIVEN BY */ /* R(J-1,J) = ( C(J) S(J) ) . */ /* ( -S(J) C(J) ) */ /* C(1) AND S(1) ARE NOT REFERENCED. */ /* N MUST BE AT LEAST 1. IF N=1 THEN AN IMMEDIATE RETURN TO */ /* THE CALLING PROGRAM IS MADE. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --x; --s; --c__; /* Function Body */ if (*n == 1) { return 0; } i__ = *n; i__1 = *n; for (ii = 2; ii <= i__1; ++ii) { im1 = i__ - 1; w = x[im1]; x[im1] = c__[i__] * w + s[i__] * x[i__]; x[i__] = c__[i__] * x[i__] - s[i__] * w; i__ = im1; /* L20: */ } return 0; } /* f01lzx_ */ /* Subroutine */ int f01lzy_(n, c__, s, x, y) integer *n; doublereal *c__, *s, *x, *y; { /* System generated locals */ integer i__1; /* Local variables */ static integer i__; static doublereal w; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (PLROT8) */ /* F01LZY FORMS THE N*2 MATRIX */ /* Z = ( X Y )*( C -S ) , */ /* ( S C ) */ /* WHERE X AND Y ARE N ELEMENT VECTORS, C=COS(THETA) AND */ /* S=SIN(THETA). */ /* THE FIRST COLUMN OF Z IS OVERWRITTEN ON X AND THE SECOND */ /* COLUMN OF Z IS OVERWRITTEN ON Y. */ /* N MUST BE AT LEAST 1. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --y; --x; /* Function Body */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { w = x[i__]; x[i__] = *c__ * w + *s * y[i__]; y[i__] = *c__ * y[i__] - *s * w; /* L20: */ } return 0; } /* f01lzy_ */ doublereal f01lzz_(a, b, small, big) doublereal *a, *b, *small, *big; { /* System generated locals */ doublereal ret_val; /* Builtin functions */ double d_sign(); /* Local variables */ static doublereal absa, absb, x; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988). */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (TANGNT) */ /* F01LZZ RETURNS THE VALUE */ /* F01LZZ = B/A . */ /* SMALL AND BIG MUST BE SUCH THAT */ /* SMALL = X02AMF AND BIG = 1.0/SMALL , */ /* WHERE X02AMF IS THE SMALL NUMBER RETURNED FROM ROUTINE */ /* X02AMF. */ /* IF B/A IS LESS THAN SMALL THEN F01LZZ IS RETURNED AS */ /* ZERO AND IF B/A IS GREATER THAN BIG THEN F01LZZ IS */ /* RETURNED AS SIGN(BIG,B). */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ ret_val = 0.; if (*b == 0.) { return ret_val; } absa = abs(*a); absb = abs(*b); x = 0.; if (absa >= 1.) { x = absa * *small; } if (absb < x) { return ret_val; } x = 0.; if (absb >= 1.) { x = absb * *small; } if (absa <= x) { goto L20; } ret_val = *b / *a; return ret_val; L20: ret_val = d_sign(big, b); return ret_val; } /* f01lzz_ */ /* Subroutine */ int f01qcf_(m, n, a, lda, zeta, ifail) integer *m, *n; doublereal *a; integer *lda; doublereal *zeta; integer *ifail; { /* Format strings */ static char fmt_99999[] = "(\002 The input parameters contained \002,\ i2,\002 error(s)\002)"; /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; /* Builtin functions */ integer s_wsfi(), do_fio(), e_wsfi(); /* Local variables */ extern /* Subroutine */ int dger_(); static integer ierr; static doublereal temp; extern integer p01abf_(); static integer k; extern /* Subroutine */ int f06frf_(), p01aby_(), dgemv_(); static integer la; static char rec[46*1]; /* Fortran I/O blocks */ static icilist io___76 = { 0, rec, 0, fmt_99999, 46, 1 }; /* MARK 14 RELEASE. NAG COPYRIGHT 1989. */ /* 1. Purpose */ /* ======= */ /* F01QCF finds the QR factorization of the real m by n, m .ge. n, */ /* matrix A, so that A is reduced to upper triangular form by means of */ /* orthogonal transformations. */ /* 2. Description */ /* =========== */ /* The m by n matrix A is factorized as */ /* A = Q*( R ) when m.gt.n, */ /* ( 0 ) */ /* A = Q*R when m = n, */ /* where Q is an m by m orthogonal matrix and R is an n by n upper */ /* triangular matrix. */ /* The factorization is obtained by Householder's method. The kth */ /* transformation matrix, Q( k ), which is used to introduce zeros into */ /* the kth column of A is given in the form */ /* Q( k ) = ( I 0 ), */ /* ( 0 T( k ) ) */ /* where */ /* T( k ) = I - u( k )*u( k )', */ /* u( k ) = ( zeta( k ) ), */ /* ( z( k ) ) */ /* zeta( k ) is a scalar and z( k ) is an ( m - k ) element vector. */ /* zeta( k ) and z( k ) are chosen to annhilate the elements below the */ /* triangular part of A. */ /* The vector u( k ) is returned in the kth element of ZETA and in the */ /* kth column of A, such that zeta( k ) is in ZETA( k ) and the elements */ /* of z( k ) are in a( k + 1, k ), ..., a( m, k ). The elements of R */ /* are returned in the upper triangular part of A. */ /* Q is given by */ /* Q = ( Q( n )*Q( n - 1 )*...*Q( 1 ) )'. */ /* 3. Parameters */ /* ========== */ /* M - INTEGER. */ /* On entry, M must specify the number of rows of A. M must be */ /* at least n. */ /* Unchanged on exit. */ /* N - INTEGER. */ /* On entry, N must specify the number of columns of A. N must */ /* be at least zero. When N = 0 then an immediate return is */ /* effected. */ /* Unchanged on exit. */ /* A - REAL array of DIMENSION ( LDA, n ). */ /* Before entry, the leading M by N part of the array A must */ /* contain the matrix to be factorized. */ /* On exit, the N by N upper triangular part of A will contain */ /* the upper triangular matrix R and the M by N strictly lower */ /* triangular part of A will contain details of the */ /* factorization as described above. */ /* LDA - INTEGER. */ /* On entry, LDA must specify the leading dimension of the */ /* array A as declared in the calling (sub) program. LDA must */ /* be at least m. */ /* Unchanged on exit. */ /* ZETA - REAL array of DIMENSION at least ( n ). */ /* On exit, ZETA( k ) contains the scalar zeta( k ) for the */ /* kth transformation. If T( k ) = I then ZETA( k ) = 0.0, */ /* otherwise ZETA( k ) contains zeta( k ) as described above */ /* and zeta( k ) is always in the range ( 1.0, sqrt( 2.0 ) ). */ /* IFAIL - INTEGER. */ /* Before entry, IFAIL must contain one of the values -1 or 0 */ /* or 1 to specify noisy soft failure or noisy hard failure or */ /* silent soft failure. ( See Chapter P01 for further details.) */ /* On successful exit IFAIL will be zero, otherwise IFAIL */ /* will be set to -1 indicating that an input parameter has */ /* been incorrectly set. See the next section for further */ /* details. */ /* 4. Diagnostic Information */ /* ====================== */ /* IFAIL = -1 */ /* One or more of the following conditions holds: */ /* M .lt. N */ /* N .lt. 0 */ /* LDA .lt. M */ /* If on entry, IFAIL was either -1 or 0 then further diagnostic */ /* information will be output on the error message channel. ( See */ /* routine X04AAF. ) */ /* 5. Further information */ /* =================== */ /* Following the use of this routine the operations */ /* B := Q*B and B := Q'*B, */ /* where B is an m by k matrix, can be performed by calls to the */ /* NAG Library routine F01QDF. The operation B := Q*B can be obtained */ /* by the call: */ /* IFAIL = 0 */ /* CALL F01QDF( 'No transpose', 'Separate', M, N, A, LDA, ZETA, */ /* $ K, B, LDB, WORK, IFAIL ) */ /* and B := Q'*B can be obtained by the call: */ /* IFAIL = 0 */ /* CALL F01QDF( 'Transpose', 'Separate', M, N, A, LDA, ZETA, */ /* $ K, B, LDB, WORK, IFAIL ) */ /* In both cases WORK must be a k element array that is used as */ /* workspace. If B is a one-dimensional array (single column) then the */ /* parameter LDB can be replaced by M. See routine F01QDF for further */ /* details. */ /* The first k columns of the orthogonal matrix Q can either be obtained */ /* by setting B to the first k columns of the unit matrix and using the */ /* first of the above two calls, or by calling the NAG Library routine */ /* F01QEF, which overwrites the k columns of Q on the first k columns of */ /* the array A. Q is obtained by the call: */ /* CALL F01QEF( 'Separate', M, N, K, A, LDA, ZETA, WORK, IFAIL ) */ /* As above WORK must be a k element array. If K is larger than N, then */ /* A must have been declared to have at least K columns. */ /* Operations involving the matrix R can readily be performed by the */ /* Level 2 BLAS routines DTRSV and DTRMV (see Chapter F06), but note */ /* that no test for near singularity of R is incorporated in DTRSV . */ /* If R is singular, or nearly singular then the NAG Library routine */ /* F02WUF can be used to determine the singular value decomposition */ /* of R. */ /* Nag Fortran 77 Auxiliary linear algebra routine. */ /* -- Written on 21-December-1985. */ /* Sven Hammarling, Nag Central Office. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Check the input parameters. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; --zeta; /* Function Body */ ierr = 0; if (*m < *n) { p01aby_(m, "M", ifail, &ierr, "F01QCF", 1L, 6L); } if (*n < 0) { p01aby_(n, "N", ifail, &ierr, "F01QCF", 1L, 6L); } if (*lda < *m) { p01aby_(lda, "LDA", ifail, &ierr, "F01QCF", 3L, 6L); } if (ierr > 0) { s_wsfi(&io___76); do_fio(&c__1, (char *)&ierr, (ftnlen)sizeof(integer)); e_wsfi(); *ifail = p01abf_(ifail, &c_n1, "F01QCF", &c__1, rec, 6L, 46L); return 0; } /* Perform the factorization. */ if (*n == 0) { *ifail = p01abf_(ifail, &c__0, "F01QCF", &c__0, rec, 6L, 46L); return 0; } la = *lda; /* Computing MIN */ i__2 = *m - 1; i__1 = min(i__2,*n); for (k = 1; k <= i__1; ++k) { /* Use a Householder reflection to zero the kth column of A. */ /* First set up the reflection. */ i__2 = *m - k; f06frf_(&i__2, &a[k + k * a_dim1], &a[k + 1 + k * a_dim1], &c__1, & c_b114, &zeta[k]); if (zeta[k] > 0. && k < *n) { if (k + 1 == *n) { la = *m - k + 1; } /* Temporarily store beta and put zeta( k ) in a( k, k ). */ temp = a[k + k * a_dim1]; a[k + k * a_dim1] = zeta[k]; /* We now perform the operation A := Q( k )*A. */ /* Let B denote the bottom ( m - k + 1 ) by ( n - k ) part */ /* of A. */ /* First form work = B'*u. ( work is stored in the e lements */ /* ZETA( k + 1 ), ..., ZETA( n ). ) */ i__2 = *m - k + 1; i__3 = *n - k; dgemv_("Transpose", &i__2, &i__3, &c_b15, &a[k + (k + 1) * a_dim1] , &la, &a[k + k * a_dim1], &c__1, &c_b114, &zeta[k + 1], & c__1, 9L); /* Now form B := B - u*work'. */ i__2 = *m - k + 1; i__3 = *n - k; dger_(&i__2, &i__3, &c_b24, &a[k + k * a_dim1], &c__1, &zeta[k + 1], &c__1, &a[k + (k + 1) * a_dim1], &la); /* Restore beta. */ a[k + k * a_dim1] = temp; } /* L20: */ } /* Set the final ZETA when m.eq.n. */ if (*m == *n) { zeta[*n] = 0.; } *ifail = p01abf_(ifail, &c__0, "F01QCF", &c__0, rec, 6L, 46L); return 0; /* End of F01QCF. ( SGEQR ) */ } /* f01qcf_ */ /* Subroutine */ int f01qdf_(trans, wheret, m, n, a, lda, zeta, ncolb, b, ldb, work, ifail, trans_len, wheret_len) char *trans, *wheret; integer *m, *n; doublereal *a; integer *lda; doublereal *zeta; integer *ncolb; doublereal *b; integer *ldb; doublereal *work; integer *ifail; ftnlen trans_len; ftnlen wheret_len; { /* Format strings */ static char fmt_99999[] = "(\002 The input parameters contained \002,\ i2,\002 error(s)\002)"; /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2; /* Builtin functions */ integer s_wsfi(), do_fio(), e_wsfi(); /* Local variables */ extern /* Subroutine */ int dger_(); static integer ierr; static doublereal temp; extern integer p01abf_(); static integer k; extern /* Subroutine */ int p01abw_(), p01aby_(), dgemv_(); static doublereal zetak; static integer lb, kk; static char rec[46*1]; /* Fortran I/O blocks */ static icilist io___82 = { 0, rec, 0, fmt_99999, 46, 1 }; /* MARK 14 RELEASE. NAG COPYRIGHT 1989. */ /* On entry, LDB must specify the leading dimension of the */ /* array B as declared in the calling (sub) program. LDB must */ /* be at least m. */ /* Unchanged on exit. */ /* WORK - REAL array of DIMENSION at least ( ncolb ). */ /* Used as internal workspace. */ /* IFAIL - INTEGER. */ /* Before entry, IFAIL must contain one of the values -1 or 0 */ /* or 1 to specify noisy soft failure or noisy hard failure or */ /* silent soft failure. ( See Chapter P01 for further details.) */ /* On successful exit IFAIL will be zero, otherwise IFAIL */ /* will be set to -1 indicating that an input parameter has */ /* been incorrectly set. See the next section for further */ /* details. */ /* 4. Diagnostic Information */ /* ====================== */ /* IFAIL = -1 */ /* One or more of the following conditions holds: */ /* TRANS .ne. 'N' or 'n' or 'T' or 't' or 'C' or 'c' */ /* WHERET .ne. 'I' or 'i' or 'S' or 's' */ /* M .lt. N */ /* N .lt. 0 */ /* LDA .lt. M */ /* NCOLB .lt. 0 */ /* LDB .lt. M */ /* If on entry, IFAIL was either -1 or 0 then further diagnostic */ /* information will be output on the error message channel. ( See */ /* routine X04AAF. ) */ /* Nag Fortran 77 Auxiliary linear algebra routine. */ /* -- Written on 13-November-1987. */ /* Sven Hammarling, Nag Central Office. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Check the input parameters. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; --zeta; b_dim1 = *ldb; b_offset = b_dim1 + 1; b -= b_offset; --work; /* Function Body */ ierr = 0; if (*(unsigned char *)trans != 'N' && *(unsigned char *)trans != 'n' && *( unsigned char *)trans != 'T' && *(unsigned char *)trans != 't' && *(unsigned char *)trans != 'C' && *(unsigned char *)trans != 'c') { p01abw_(trans, "TRANS", ifail, &ierr, "F01QDF", 1L, 5L, 6L); } if (*(unsigned char *)wheret != 'I' && *(unsigned char *)wheret != 'i' && *(unsigned char *)wheret != 'S' && *(unsigned char *)wheret != 's') { p01abw_(wheret, "WHERET", ifail, &ierr, "F01QDF", 1L, 6L, 6L); } if (*m < *n) { p01aby_(m, "M", ifail, &ierr, "F01QDF", 1L, 6L); } if (*n < 0) { p01aby_(n, "N", ifail, &ierr, "F01QDF", 1L, 6L); } if (*lda < *m) { p01aby_(lda, "LDA", ifail, &ierr, "F01QDF", 3L, 6L); } if (*ncolb < 0) { p01aby_(ncolb, "NCOLB", ifail, &ierr, "F01QDF", 5L, 6L); } if (*ldb < *m) { p01aby_(ldb, "LDB", ifail, &ierr, "F01QDF", 3L, 6L); } if (ierr > 0) { s_wsfi(&io___82); do_fio(&c__1, (char *)&ierr, (ftnlen)sizeof(integer)); e_wsfi(); *ifail = p01abf_(ifail, &c_n1, "F01QDF", &c__1, rec, 6L, 46L); return 0; } /* Perform the transformation. */ if (min(*n,*ncolb) == 0) { *ifail = p01abf_(ifail, &c__0, "F01QDF", &c__0, rec, 6L, 46L); return 0; } lb = *ldb; i__1 = *n; for (kk = 1; kk <= i__1; ++kk) { if (*(unsigned char *)trans == 'T' || *(unsigned char *)trans == 't' || *(unsigned char *)trans == 'C' || *(unsigned char *)trans == 'c') { /* Q'*B = Q( n )*...*Q( 2 )*Q( 1 )*B, */ k = kk; } else { /* Q*B = Q( 1 )'*Q( 2 )'*...*Q( n )'*B, */ k = *n + 1 - kk; } if (*(unsigned char *)wheret == 'S' || *(unsigned char *)wheret == 's') { zetak = zeta[k]; } else { zetak = a[k + k * a_dim1]; } if (zetak > 0.) { temp = a[k + k * a_dim1]; a[k + k * a_dim1] = zetak; if (*ncolb == 1) { lb = *m - k + 1; } /* Let C denote the bottom ( m - k + 1 ) by ncolb part o f B. */ /* First form work = C'*u. */ i__2 = *m - k + 1; dgemv_("Transpose", &i__2, ncolb, &c_b15, &b[k + b_dim1], &lb, &a[ k + k * a_dim1], &c__1, &c_b114, &work[1], &c__1, 9L); /* Now form C := C - u*work'. */ i__2 = *m - k + 1; dger_(&i__2, ncolb, &c_b24, &a[k + k * a_dim1], &c__1, &work[1], & c__1, &b[k + b_dim1], &lb); /* Restore the diagonal element of A. */ a[k + k * a_dim1] = temp; } /* L20: */ } *ifail = p01abf_(ifail, &c__0, "F01QDF", &c__0, rec, 6L, 46L); return 0; /* End of F01QDF. ( SGEAPQ ) */ } /* f01qdf_ */ /* Subroutine */ int f02agf_(a, ia, n, rr, ri, vr, ivr, vi, ivi, intger, ifail) doublereal *a; integer *ia, *n; doublereal *rr, *ri, *vr; integer *ivr; doublereal *vi; integer *ivi, *intger, *ifail; { /* System generated locals */ integer a_dim1, a_offset, vi_dim1, vi_offset, vr_dim1, vr_offset, i__1, i__2; doublereal d__1, d__2; /* Builtin functions */ double sqrt(); /* Local variables */ static doublereal term, c__, d__; extern /* Subroutine */ int f01akf_(); static integer i__, j, k, l; extern integer p01abf_(); extern /* Subroutine */ int f01apf_(), f02aqf_(), f01atf_(), f01auf_(); extern integer x02bhf_(); static char p01rec[1*1]; extern doublereal x02ajf_(); static integer isave, ib; static doublereal machep, max__, sum; /* MARK 13 RE-ISSUE. NAG COPYRIGHT 1988. */ /* MARK 14A REVISED. IER-685 (DEC 1989). */ /* EIGENVALUES AND EIGENVECTORS OF REAL UNSYMMETRIC MATRIX */ /* 1ST AUGUST 1971 */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --intger; --ri; --rr; a_dim1 = *ia; a_offset = a_dim1 + 1; a -= a_offset; vr_dim1 = *ivr; vr_offset = vr_dim1 + 1; vr -= vr_offset; vi_dim1 = *ivi; vi_offset = vi_dim1 + 1; vi -= vi_offset; /* Function Body */ isave = *ifail; *ifail = 1; machep = x02ajf_(); ib = x02bhf_(); f01atf_(n, &ib, &a[a_offset], ia, &k, &l, &rr[1]); f01akf_(n, &k, &l, &a[a_offset], ia, &intger[1]); f01apf_(n, &k, &l, &intger[1], &a[a_offset], ia, &vr[vr_offset], ivr); f01auf_(n, &k, &l, n, &rr[1], &vr[vr_offset], ivr); f02aqf_(n, &c__1, n, &machep, &a[a_offset], ia, &vr[vr_offset], ivr, &rr[ 1], &ri[1], &intger[1], ifail); if (*ifail == 0) { goto L20; } *ifail = p01abf_(&isave, ifail, "F02AGF", &c__0, p01rec, 6L, 1L); return 0; L20: i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (ri[i__] == 0.) { goto L60; } if (ri[i__] > 0.) { goto L100; } i__2 = *n; for (j = 1; j <= i__2; ++j) { vr[j + i__ * vr_dim1] = vr[j + (i__ - 1) * vr_dim1]; vi[j + i__ * vi_dim1] = -vi[j + (i__ - 1) * vi_dim1]; /* L40: */ } goto L140; L60: i__2 = *n; for (j = 1; j <= i__2; ++j) { vi[j + i__ * vi_dim1] = 0.; /* L80: */ } goto L140; L100: i__2 = *n; for (j = 1; j <= i__2; ++j) { vi[j + i__ * vi_dim1] = vr[j + (i__ + 1) * vr_dim1]; /* L120: */ } L140: ; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { sum = 0.; max__ = 0.; i__2 = *n; for (j = 1; j <= i__2; ++j) { if ((d__1 = vr[j + i__ * vr_dim1], abs(d__1)) <= max__) { goto L160; } max__ = (d__1 = vr[j + i__ * vr_dim1], abs(d__1)); L160: if ((d__1 = vi[j + i__ * vi_dim1], abs(d__1)) <= max__) { goto L180; } max__ = (d__1 = vi[j + i__ * vi_dim1], abs(d__1)); L180: ; } i__2 = *n; for (j = 1; j <= i__2; ++j) { vr[j + i__ * vr_dim1] /= max__; vi[j + i__ * vi_dim1] /= max__; /* L200: */ } max__ = 0.; i__2 = *n; for (j = 1; j <= i__2; ++j) { /* Computing 2nd power */ d__1 = vr[j + i__ * vr_dim1]; /* Computing 2nd power */ d__2 = vi[j + i__ * vi_dim1]; term = d__1 * d__1 + d__2 * d__2; sum += term; if (term <= max__) { goto L220; } max__ = term; c__ = vr[j + i__ * vr_dim1]; d__ = -vi[j + i__ * vi_dim1]; L220: /* L240: */ ; } /* Computing 2nd power */ d__1 = c__; /* Computing 2nd power */ d__2 = d__; sum *= d__1 * d__1 + d__2 * d__2; sum = sqrt(sum); i__2 = *n; for (j = 1; j <= i__2; ++j) { term = vr[j + i__ * vr_dim1]; vr[j + i__ * vr_dim1] = (vr[j + i__ * vr_dim1] * c__ - vi[j + i__ * vi_dim1] * d__) / sum; vi[j + i__ * vi_dim1] = (d__ * term + c__ * vi[j + i__ * vi_dim1]) / sum; /* L260: */ } /* L280: */ } return 0; } /* f02agf_ */ /* Subroutine */ int f02aqf_(n, low, upp, machep, h__, ih, vecs, ivecs, wr, wi, cnt, ifail) integer *n, *low, *upp; doublereal *machep, *h__; integer *ih; doublereal *vecs; integer *ivecs; doublereal *wr, *wi; integer *cnt, *ifail; { /* System generated locals */ integer h_dim1, h_offset, vecs_dim1, vecs_offset, i__1, i__2, i__3; doublereal d__1, d__2, d__3; /* Builtin functions */ double sqrt(); /* Local variables */ static doublereal norm; extern doublereal a02abf_(); extern /* Subroutine */ int a02acf_(); static integer i__, j, k, l, m; static doublereal p, q, r__, s, t, u, w, x, y, z__; static char p01rec[1*1]; extern doublereal x02ajf_(), x02akf_(); extern integer p01abf_(); extern /* Subroutine */ int dgemv_(); static integer isave, i1, m2, m3, na, ii; static doublereal ra, sa; static integer en, jj, kk, ll, mm; static doublereal vi, vr; static integer na1; static logical notlas; static integer en2, nhs, itn, its, low1, upp1; /* MARK 13 RE-ISSUE. NAG COPYRIGHT 1988. */ /* HQR2 */ /* FINDS THE EIGENVALUES AND EIGENVECTORS OF A REAL MATRIX */ /* WHICH HAS BEEN REDUCED TO UPPER HESSENBERG FORM IN THE ARRAY */ /* H(N,N) WITH THE ACCUMULATED TRANSFORMATIONS STORED IN */ /* THE ARRAY VECS(N,N). THE REAL AND IMAGINARY PARTS OF THE */ /* EIGENVALUES ARE FORMED IN THE ARRAYS WR, WI(N) AND THE */ /* EIGENVECTORS ARE FORMED IN THE ARRAY VECS(N,N) WHERE */ /* ONLY ONE COMPLEX VECTOR, CORRESPONDING TO THE ROOT WITH */ /* POSITIVE IMAGINARY PART, IS FORMED FOR A COMPLEX PAIR. LOW */ /* AND UPP ARE TWO INTEGERS PRODUCED IN BALANCING WHERE */ /* EIGENVALUES ARE ISOLATED IN POSITIONS 1 TO LOW-1 AND UPP+1 */ /* TO N. IF BALANCING IS NOT USED LOW=1, UPP=N. MACHEP IS THE */ /* RELATIVE MACHINE PRECISION. THE SUBROUTINE WILL FAIL IF */ /* ALL EIGENVALUES TAKE MORE THAN 30*N ITERATIONS. */ /* 1ST DECEMBER 1971 */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --cnt; --wi; --wr; h_dim1 = *ih; h_offset = h_dim1 + 1; h__ -= h_offset; vecs_dim1 = *ivecs; vecs_offset = vecs_dim1 + 1; vecs -= vecs_offset; /* Function Body */ isave = *ifail; /* COMPUTE MATRIX NORM */ norm = 0.; k = 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *n; for (j = k; j <= i__2; ++j) { norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1)); /* L20: */ } k = i__; /* L40: */ } nhs = *n * (*n + 1) / 2 + *n - 1; /* ISOLATED ROOTS */ if (*low <= 1) { goto L80; } j = *low - 1; i__1 = j; for (i__ = 1; i__ <= i__1; ++i__) { wr[i__] = h__[i__ + i__ * h_dim1]; wi[i__] = 0.; cnt[i__] = 0; /* L60: */ } L80: if (*upp >= *n) { goto L120; } j = *upp + 1; i__1 = *n; for (i__ = j; i__ <= i__1; ++i__) { wr[i__] = h__[i__ + i__ * h_dim1]; wi[i__] = 0.; cnt[i__] = 0; /* L100: */ } L120: en = *upp; t = 0.; itn = *n * 30; L140: if (en < *low) { goto L880; } its = 0; na = en - 1; /* LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */ L160: if (*low + 1 > en) { goto L200; } low1 = *low + 1; i__1 = en; for (ll = low1; ll <= i__1; ++ll) { l = en + low1 - ll; s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l + l * h_dim1], abs(d__2)); if (s < x02akf_() / x02ajf_()) { s = norm / (doublereal) nhs; } if ((d__1 = h__[l + (l - 1) * h_dim1], abs(d__1)) <= *machep * s) { goto L220; } /* L180: */ } L200: l = *low; L220: x = h__[en + en * h_dim1]; if (l == en) { goto L740; } y = h__[na + na * h_dim1]; w = h__[en + na * h_dim1] * h__[na + en * h_dim1]; if (l == na) { goto L760; } if (itn <= 0) { goto L1500; } /* FORM SHIFT */ if (its != 10 && its != 20) { goto L280; } t += x; if (*low > en) { goto L260; } i__1 = en; for (i__ = *low; i__ <= i__1; ++i__) { h__[i__ + i__ * h_dim1] -= x; /* L240: */ } L260: s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + (en - 2) * h_dim1], abs(d__2)); x = s * .75; y = x; /* Computing 2nd power */ d__1 = s; w = d__1 * d__1 * -.4375; L280: ++its; --itn; /* LOOK FOR TWO CONSECUTIVE SMALL SUB-DIAGONAL ELEMENTS */ if (l > en - 2) { goto L320; } en2 = en - 2; i__1 = en2; for (mm = l; mm <= i__1; ++mm) { m = l + en2 - mm; z__ = h__[m + m * h_dim1]; r__ = x - z__; s = y - z__; p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) * h_dim1]; q = h__[m + 1 + (m + 1) * h_dim1] - z__ - r__ - s; r__ = h__[m + 2 + (m + 1) * h_dim1]; s = abs(p) + abs(q) + abs(r__); p /= s; q /= s; r__ /= s; if (m == l) { goto L320; } if ((d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q) + abs(r__) ) <= *machep * abs(p) * ((d__2 = h__[m - 1 + (m - 1) * h_dim1] , abs(d__2)) + abs(z__) + (d__3 = h__[m + 1 + (m + 1) * h_dim1], abs(d__3)))) { goto L320; } /* L300: */ } L320: m2 = m + 2; if (m2 > en) { goto L360; } i__1 = en; for (i__ = m2; i__ <= i__1; ++i__) { h__[i__ + (i__ - 2) * h_dim1] = 0.; /* L340: */ } L360: m3 = m + 3; if (m3 > en) { goto L400; } i__1 = en; for (i__ = m3; i__ <= i__1; ++i__) { h__[i__ + (i__ - 3) * h_dim1] = 0.; /* L380: */ } L400: if (m > na) { goto L720; } i__1 = na; for (k = m; k <= i__1; ++k) { notlas = k != na; if (k == m) { goto L420; } p = h__[k + (k - 1) * h_dim1]; q = h__[k + 1 + (k - 1) * h_dim1]; r__ = 0.; if (notlas) { r__ = h__[k + 2 + (k - 1) * h_dim1]; } x = abs(p) + abs(q) + abs(r__); if (x == 0.) { goto L700; } p /= x; q /= x; r__ /= x; L420: /* Computing 2nd power */ d__1 = p; /* Computing 2nd power */ d__2 = q; /* Computing 2nd power */ d__3 = r__; s = sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3); if (p < 0.) { s = -s; } if (k != m) { goto L440; } if (l != m) { h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1]; } goto L460; L440: h__[k + (k - 1) * h_dim1] = -s * x; L460: p += s; x = p / s; y = q / s; z__ = r__ / s; q /= p; r__ /= p; /* ROW MODIFICATION */ if (notlas) { goto L500; } i__2 = *n; for (j = k; j <= i__2; ++j) { p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1]; h__[k + 1 + j * h_dim1] -= p * y; h__[k + j * h_dim1] -= p * x; /* L480: */ } goto L540; L500: i__2 = *n; for (j = k; j <= i__2; ++j) { p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[ k + 2 + j * h_dim1]; h__[k + 2 + j * h_dim1] -= p * z__; h__[k + 1 + j * h_dim1] -= p * y; h__[k + j * h_dim1] -= p * x; /* L520: */ } L540: j = en; if (k + 3 < en) { j = k + 3; } /* COLUMN MODIFICATION */ if (notlas) { goto L580; } i__2 = j; for (i__ = 1; i__ <= i__2; ++i__) { p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1]; h__[i__ + (k + 1) * h_dim1] -= p * q; h__[i__ + k * h_dim1] -= p; /* L560: */ } goto L620; L580: i__2 = j; for (i__ = 1; i__ <= i__2; ++i__) { p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] + z__ * h__[i__ + (k + 2) * h_dim1]; h__[i__ + (k + 2) * h_dim1] -= p * r__; h__[i__ + (k + 1) * h_dim1] -= p * q; h__[i__ + k * h_dim1] -= p; /* L600: */ } /* ACCUMULATE TRANSFORMATIONS */ L620: if (*low > *upp) { goto L700; } if (notlas) { goto L660; } i__2 = *upp; for (i__ = *low; i__ <= i__2; ++i__) { p = x * vecs[i__ + k * vecs_dim1] + y * vecs[i__ + (k + 1) * vecs_dim1]; vecs[i__ + (k + 1) * vecs_dim1] -= p * q; vecs[i__ + k * vecs_dim1] -= p; /* L640: */ } goto L700; L660: i__2 = *upp; for (i__ = *low; i__ <= i__2; ++i__) { p = x * vecs[i__ + k * vecs_dim1] + y * vecs[i__ + (k + 1) * vecs_dim1] + z__ * vecs[i__ + (k + 2) * vecs_dim1]; vecs[i__ + (k + 2) * vecs_dim1] -= p * r__; vecs[i__ + (k + 1) * vecs_dim1] -= p * q; vecs[i__ + k * vecs_dim1] -= p; /* L680: */ } L700: ; } L720: goto L160; /* ONE ROOT FOUND */ L740: wr[en] = x + t; h__[en + en * h_dim1] = wr[en]; wi[en] = 0.; cnt[en] = its; en = na; goto L140; /* TWO ROOTS FOUND */ L760: p = (y - x) / 2.; /* Computing 2nd power */ d__1 = p; q = d__1 * d__1 + w; z__ = sqrt((abs(q))); x += t; h__[en + en * h_dim1] = x; h__[na + na * h_dim1] = y + t; cnt[en] = -its; cnt[na] = its; if (q < 0.) { goto L840; } /* REAL PAIR */ if (p < 0.) { z__ = p - z__; } if (p > 0.) { z__ = p + z__; } wr[na] = x + z__; wr[en] = wr[na]; if (z__ != 0.) { wr[en] = x - w / z__; } wi[na] = 0.; wi[en] = 0.; x = h__[en + na * h_dim1]; r__ = a02abf_(&x, &z__); p = x / r__; q = z__ / r__; /* ROW MODIFICATION */ i__1 = *n; for (j = na; j <= i__1; ++j) { z__ = h__[na + j * h_dim1]; h__[na + j * h_dim1] = q * z__ + p * h__[en + j * h_dim1]; h__[en + j * h_dim1] = q * h__[en + j * h_dim1] - p * z__; /* L780: */ } /* COLUMN MODIFICATION */ i__1 = en; for (i__ = 1; i__ <= i__1; ++i__) { z__ = h__[i__ + na * h_dim1]; h__[i__ + na * h_dim1] = q * z__ + p * h__[i__ + en * h_dim1]; h__[i__ + en * h_dim1] = q * h__[i__ + en * h_dim1] - p * z__; /* L800: */ } /* ACCUMULATE TRANSFORMATIONS */ i__1 = *upp; for (i__ = *low; i__ <= i__1; ++i__) { z__ = vecs[i__ + na * vecs_dim1]; vecs[i__ + na * vecs_dim1] = q * z__ + p * vecs[i__ + en * vecs_dim1]; vecs[i__ + en * vecs_dim1] = q * vecs[i__ + en * vecs_dim1] - p * z__; /* L820: */ } goto L860; /* COMPLEX PAIR */ L840: wr[na] = x + p; wr[en] = x + p; wi[na] = z__; wi[en] = -z__; L860: en += -2; goto L140; /* ALL ROOTS FOUND NOW BACKSUBSTITUTE */ L880: if (norm == 0.) { goto L1480; } norm *= *machep; /* BACKSUBSTITUTION */ i__1 = *n; for (kk = 1; kk <= i__1; ++kk) { en = *n + 1 - kk; p = wr[en]; q = wi[en]; na = en - 1; if (q != 0.) { goto L1120; } /* REAL VECTOR */ h__[en + en * h_dim1] = 1.; if (na < 1) { goto L1340; } i__2 = na; for (ii = 1; ii <= i__2; ++ii) { i__ = na + 1 - ii; i1 = i__ - 1; w = h__[i__ + i__ * h_dim1] - p; r__ = h__[i__ + en * h_dim1]; if (wi[i__] >= 0.) { goto L900; } z__ = w; s = r__; goto L1100; L900: if (wi[i__] > 0.) { goto L1020; } /* MODIFICATION TO STOP OVERFLOW */ if (w != 0.) { goto L940; } if (abs(r__) < norm * 10.) { goto L960; } r__ = -r__; i__3 = en; for (j = 1; j <= i__3; ++j) { h__[j + en * h_dim1] *= norm; /* L920: */ } goto L980; L940: r__ = -r__ / w; goto L980; L960: r__ = -r__ / norm; L980: h__[i__ + en * h_dim1] = r__; if (i1 == 0) { goto L1100; } i__3 = i1; for (j = 1; j <= i__3; ++j) { h__[j + en * h_dim1] += h__[j + i__ * h_dim1] * r__; /* L1000: */ } goto L1100; /* SOLVE REAL EQUATIONS */ L1020: x = h__[i__ + (i__ + 1) * h_dim1]; y = h__[i__ + 1 + i__ * h_dim1]; /* Computing 2nd power */ d__1 = wr[i__] - p; /* Computing 2nd power */ d__2 = wi[i__]; q = d__1 * d__1 + d__2 * d__2; t = (x * s - z__ * r__) / q; h__[i__ + en * h_dim1] = t; if (abs(x) > abs(z__)) { goto L1040; } r__ = (-s - y * t) / z__; goto L1060; L1040: r__ = (-r__ - w * t) / x; L1060: h__[i__ + 1 + en * h_dim1] = r__; if (i1 == 0) { goto L1100; } i__3 = i1; for (j = 1; j <= i__3; ++j) { h__[j + en * h_dim1] = h__[j + en * h_dim1] + h__[j + (i__ + 1) * h_dim1] * r__ + h__[j + i__ * h_dim1] * t; /* L1080: */ } L1100: ; } /* END REAL VECTOR */ goto L1340; L1120: if (q > 0.) { goto L1340; } /* COMPLEX VECTOR ASSOCIATED WITH LAMBDA=P-I*Q */ if ((d__1 = h__[en + na * h_dim1], abs(d__1)) <= (d__2 = h__[na + en * h_dim1], abs(d__2))) { goto L1140; } r__ = q / h__[en + na * h_dim1]; s = -(h__[en + en * h_dim1] - p) / h__[en + na * h_dim1]; goto L1160; L1140: d__1 = -h__[na + en * h_dim1]; d__2 = h__[na + na * h_dim1] - p; a02acf_(&c_b114, &d__1, &d__2, &q, &r__, &s); L1160: h__[en + na * h_dim1] = 0.; h__[en + en * h_dim1] = 1.; h__[na + na * h_dim1] = r__; h__[na + en * h_dim1] = s; if (na < 2) { goto L1340; } na1 = na - 1; i__2 = na1; for (j = 1; j <= i__2; ++j) { h__[j + en * h_dim1] += h__[j + na * h_dim1] * s; h__[j + na * h_dim1] *= r__; /* L1180: */ } i__2 = na1; for (ii = 1; ii <= i__2; ++ii) { i__ = na1 + 1 - ii; i1 = i__ - 1; w = h__[i__ + i__ * h_dim1] - p; ra = h__[i__ + na * h_dim1]; sa = h__[i__ + en * h_dim1]; if (wi[i__] >= 0.) { goto L1200; } z__ = w; r__ = ra; s = sa; goto L1320; L1200: if (wi[i__] == 0.) { goto L1280; } /* SOLVE COMPLEX EQUATIONS */ x = h__[i__ + (i__ + 1) * h_dim1]; y = h__[i__ + 1 + i__ * h_dim1]; /* Computing 2nd power */ d__1 = wr[i__] - p; /* Computing 2nd power */ d__2 = wi[i__]; /* Computing 2nd power */ d__3 = q; vr = d__1 * d__1 + d__2 * d__2 - d__3 * d__3; vi = (wr[i__] - p) * 2. * q; if (vr == 0. && vi == 0.) { vr = *machep * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(z__)); } d__1 = x * r__ - z__ * ra + q * sa; d__2 = x * s - z__ * sa - q * ra; a02acf_(&d__1, &d__2, &vr, &vi, &t, &u); if (abs(x) <= abs(z__) + abs(q)) { goto L1220; } r__ = (-ra - w * t + q * u) / x; s = (-sa - w * u - q * t) / x; goto L1240; L1220: d__1 = -r__ - y * t; d__2 = -s - y * u; a02acf_(&d__1, &d__2, &z__, &q, &r__, &s); L1240: h__[i__ + na * h_dim1] = t; h__[i__ + en * h_dim1] = u; h__[i__ + 1 + na * h_dim1] = r__; h__[i__ + 1 + en * h_dim1] = s; if (i1 == 0) { goto L1320; } i__3 = i1; for (j = 1; j <= i__3; ++j) { h__[j + na * h_dim1] = h__[j + na * h_dim1] + h__[j + (i__ + 1) * h_dim1] * r__ + h__[j + i__ * h_dim1] * t; h__[j + en * h_dim1] = h__[j + en * h_dim1] + h__[j + (i__ + 1) * h_dim1] * s + h__[j + i__ * h_dim1] * u; /* L1260: */ } goto L1320; L1280: d__1 = -ra; d__2 = -sa; a02acf_(&d__1, &d__2, &w, &q, &r__, &s); h__[i__ + na * h_dim1] = r__; h__[i__ + en * h_dim1] = s; if (i1 == 0) { goto L1320; } i__3 = i1; for (j = 1; j <= i__3; ++j) { h__[j + na * h_dim1] += h__[j + i__ * h_dim1] * r__; h__[j + en * h_dim1] += h__[j + i__ * h_dim1] * s; /* L1300: */ } L1320: ; } /* END COMPLEX VECTOR */ L1340: ; } /* END BACKSUBSTITUTION */ /* VECTORS OF ISOLATED ROOTS */ low1 = *low - 1; upp1 = *upp + 1; i__1 = *n; for (j = 1; j <= i__1; ++j) { m = min(j,low1); if (m < 1) { goto L1380; } i__2 = m; for (i__ = 1; i__ <= i__2; ++i__) { vecs[i__ + j * vecs_dim1] = h__[i__ + j * h_dim1]; /* L1360: */ } L1380: if (upp1 > j) { goto L1420; } i__2 = j; for (i__ = upp1; i__ <= i__2; ++i__) { vecs[i__ + j * vecs_dim1] = h__[i__ + j * h_dim1]; /* L1400: */ } L1420: ; } /* MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */ /* VECTORS OF ORIGINAL FULL MATRIX */ i__1 = *n; for (jj = *low; jj <= i__1; ++jj) { j = *low + *n - jj; m = min(j,*upp); i__2 = *upp; for (i__ = *low; i__ <= i__2; ++i__) { vecs[i__ + j * vecs_dim1] = vecs[i__ + m * vecs_dim1] * h__[m + j * h_dim1]; /* L1440: */ } --m; if (m + 1 >= *low) { i__2 = *upp - *low + 1; i__3 = m - *low + 1; dgemv_("N", &i__2, &i__3, &c_b15, &vecs[*low + *low * vecs_dim1], ivecs, &h__[*low + j * h_dim1], &c__1, &c_b15, &vecs[*low + j * vecs_dim1], &c__1, 1L); } /* L1460: */ } L1480: *ifail = 0; return 0; L1500: *ifail = p01abf_(&isave, &c__1, "F02AQF", &c__0, p01rec, 6L, 1L); return 0; } /* f02aqf_ */ /* Subroutine */ int f02szf_(n, d__, e, sv, wantb, b, wanty, y, nry, ly, wantz, z__, nrz, ncz, work1, work2, work3, ifail) integer *n; doublereal *d__, *e, *sv; logical *wantb; doublereal *b; logical *wanty; doublereal *y; integer *nry, *ly; logical *wantz; doublereal *z__; integer *nrz, *ncz; doublereal *work1, *work2, *work3; integer *ifail; { /* System generated locals */ integer y_dim1, y_offset, z_dim1, z_offset, i__1, i__2, i__3; doublereal d__1, d__2, d__3, d__4; /* Builtin functions */ double sqrt(); /* Local variables */ static integer ierr, iter; static doublereal c__, f, g; static integer i__, j, k, l; extern integer p01abf_(); static doublereal s, t, x; static char p01rec[1*1]; extern doublereal x02ajf_(), x02amf_(); static doublereal small, anorm; static integer maxit; extern /* Subroutine */ int f01lzw_(), f01lzy_(); extern doublereal f01lzz_(); static doublereal shuft; extern /* Subroutine */ int f02szz_(); static doublereal dk, dl, ek; static integer jj, kk, ll, lm1, lp1; static doublereal sqteps, rsqtps, big, eps, svi, dkm1, ekm1; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 9 REVISED. IER-328 (SEP 1981). */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 12 REVISED. IER-518 (AUG 1986). */ /* MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988). */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (SVDBID) */ /* F02SZF RETURNS PART OR ALL OF THE SINGULAR VALUE */ /* DECOMPOSITION OF THE N*N UPPER BIDIAGONAL MATRIX A. THAT */ /* IS, A IS FACTORIZED AS */ /* A = Q*DIAG(SV)*(P**T) , */ /* WHERE Q AND P ARE N*N ORTHOGONAL MATRICES AND DIAG(SV) */ /* IS AN N*N DIAGONAL MATRIX WITH NON-NEGATIVE DIAGONAL */ /* ELEMENTS SV(1),SV(2),..., SV(N), THESE BEING THE */ /* SINGULAR VALUES OF A. */ /* IF WANTB IS .TRUE. THEN B RETURNS (Q**T)*B. */ /* IF WANTY IS .TRUE. THEN Y RETURNS Y*Q. */ /* IF WANTZ IS .TRUE. THEN Z RETURNS (P**T)*Z. */ /* INPUT PARAMETERS. */ /* N - THE ORDER OF THE MATRIX. MUST BE AT LEAST 1. */ /* D - N ELEMENT VECTOR SUCH THAT D(I)=A(I,I), I=1,2,...,N. */ /* D IS UNALTERED UNLESS ROUTINE IS CALLED WITH SV=D. */ /* E - N ELEMENT VECTOR SUCH THAT E(I)=A(I-1,I), I=2,3,...,N. */ /* E(1) IS NOT REFERENCED. */ /* E IS UNALTERED UNLESS ROUTINE IS CALLED WITH WORK1=E. */ /* WANTB - MUST BE .TRUE. IF (Q**T)*B IS REQUIRED. */ /* IF WANTB IS .FALSE. THEN B IS NOT REFERENCED. */ /* B - AN N ELEMENT REAL VECTOR. */ /* WANTY - MUST BE .TRUE. IF Y*Q IS REQUIRED. */ /* IF WANTY IS .FALSE. THEN Y IS NOT REFERENCED. */ /* Y - AN LY*N REAL MATRIX. */ /* NRY - IF WANTY IS .TRUE. THEN NRY MUST BE THE ROW */ /* DIMENSION OF Y AS DECLARED IN THE CALLING */ /* PROGRAM AND MUST BE AT LEAST LY. */ /* LY - IF WANTY IS .TRUE. THEN LY MUST BE THE NUMBER */ /* OF ROWS OF Y AND MUST BE AT LEAST 1. */ /* WANTZ - MUST BE .TRUE. IF (P**T)*Z IS REQUIRED. */ /* IF WANTZ IS .FALSE. THEN Z IS NOT REFERENCED. */ /* Z - AN N*NCZ REAL MATRIX. */ /* NRZ - IF WANTZ IS .TRUE. THEN NRZ MUST BE THE ROW */ /* DIMENSION OF Z AS DECLARED IN THE CALLING */ /* PROGRAM AND MUST BE AT LEAST N. */ /* NCZ - IF WANTZ IS .TRUE. THEN NCZ MUST BE THE */ /* NUMBER OF COLUMNS OF Z AND MUST BE AT LEAST */ /* 1. */ /* IFAIL - THE USUAL FAILURE PARAMETER. IF IN DOUBT SET */ /* IFAIL TO ZERO BEFORE CALLING F02SZF. */ /* OUTPUT PARAMETERS. */ /* SV - N ELEMENT VECTOR CONTAINING THE SINGULAR */ /* VALUES OF A. THEY ARE ORDERED SO THAT */ /* SV(1).GE.SV(2).GE. ... .GE.SV(N). THE ROUTINE */ /* MAY BE CALLED WITH SV=D. */ /* B - IF WANTB IS .TRUE. THEN B WILL RETURN THE N */ /* ELEMENT VECTOR (Q**T)*B. */ /* Y - IF WANTY IS .TRUE. THEN Y WILL RETURN THE */ /* LY*N MATRIX Y*Q. */ /* Z - IF WANTZ IS .TRUE. THEN Z WILL RETURN THE N*NCZ MATRIX */ /* (P**T)*Z. */ /* IFAIL - ON NORMAL RETURN IFAIL WILL BE ZERO. */ /* IN THE UNLIKELY EVENT THAT THE QR-ALGORITHM */ /* FAILS TO FIND THE SINGULAR VALUES IN 50*N */ /* ITERATIONS THEN IFAIL WILL BE 2 OR MORE AND */ /* SUCH THAT SV(1),SV(2),..,SV(IFAIL-1) MAY NOT */ /* HAVE BEEN FOUND. SEE WORK1 BELOW. THIS */ /* FAILURE IS NOT LIKELY TO OCCUR. */ /* IF AN INPUT PARAMETER IS INCORRECTLY SUPPLIED */ /* THEN IFAIL IS SET TO UNITY. */ /* WORKSPACE PARAMETERS. */ /* WORK1 - AN N ELEMENT VECTOR. IF E IS NOT REQUIRED ON */ /* RETURN THEN THE ROUTINE MAY BE CALLED WITH */ /* WORK1=E. WORK1(1) RETURNS THE TOTAL NUMBER OF */ /* ITERATIONS TAKEN BY THE QR-ALGORITHM. IF */ /* IFAIL IS POSITIVE ON RETURN THEN THE MATRIX A */ /* IS GIVEN BY A=Q*C*(P**T) , WHERE C IS THE */ /* UPPER BIDIAGONAL MATRIX WITH SV AS ITS */ /* DIAGONAL AND WORK1 AS ITS SUPER-DIAGONAL. */ /* WORK2 */ /* WORK3 - N ELEMENT VECTORS. IF WANTZ IS .FALSE. THEN WORK2 AND */ /* WORK3 ARE NOT REFERENCED. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --work3; --work2; --work1; --b; --sv; --e; --d__; y_dim1 = *nry; y_offset = y_dim1 + 1; y -= y_offset; z_dim1 = *nrz; z_offset = z_dim1 + 1; z__ -= z_offset; /* Function Body */ ierr = *ifail; if (ierr == 0) { *ifail = 1; } if (*n < 1) { goto L500; } if (*wanty && (*nry < *ly || *ly < 1)) { goto L500; } if (*wantz && (*nrz < *n || *ncz < 1)) { goto L500; } small = x02amf_(); big = 1. / small; eps = x02ajf_(); sqteps = sqrt(eps); rsqtps = 1. / sqteps; iter = 0; k = *n; sv[1] = d__[1]; anorm = abs(d__[1]); if (*n == 1) { goto L280; } i__1 = *n; for (i__ = 2; i__ <= i__1; ++i__) { sv[i__] = d__[i__]; work1[i__] = e[i__]; /* Computing MAX */ d__3 = anorm, d__4 = (d__1 = d__[i__], abs(d__1)), d__3 = max(d__3, d__4), d__4 = (d__2 = e[i__], abs(d__2)); anorm = max(d__3,d__4); /* L20: */ } maxit = *n * 50; eps *= anorm; /* MAXIT IS THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. */ /* EPS WILL BE USED TO TEST FOR NEGLIGIBLE ELEMENTS. */ /* START MAIN LOOP. ONE SINGULAR VALUE IS FOUND FOR EACH */ /* VALUE OF K. K GOES IN OPPOSITE DIRECTION TO KK. */ i__1 = *n; for (kk = 2; kk <= i__1; ++kk) { /* NOW TEST FOR SPLITTING. L GOES IN OPPOSITE DIRECTION TO LL. */ L40: l = k; i__2 = k; for (ll = 2; ll <= i__2; ++ll) { if ((d__1 = work1[l], abs(d__1)) <= eps) { goto L240; } --l; if ((d__1 = sv[l], abs(d__1)) < eps) { goto L180; } /* L60: */ } L80: if (iter == maxit) { goto L280; } /* MAXIT QR-STEPS WITHOUT CONVERGENCE. FAILURE. */ ++iter; /* NOW DETERMINE SHIFT. */ lp1 = l + 1; dl = sv[l]; dkm1 = sv[k - 1]; dk = sv[k]; ekm1 = 0.; if (k != 2) { ekm1 = work1[k - 1]; } ek = work1[k]; f = (dkm1 - dk) * (dkm1 + dk) + (ekm1 - ek) * (ekm1 + ek); f /= ek * 2. * dkm1; g = abs(f); if (g <= rsqtps) { /* Computing 2nd power */ d__1 = f; g = sqrt(d__1 * d__1 + 1.); } if (f < 0.) { g = -g; } shuft = ek * (ek - dkm1 / (f + g)); f = (dl - dk) * (dl + dk) - shuft; x = dl * work1[lp1]; /* NOW PERFORM THE QR-STEP AND CHASE ZEROS. */ i__2 = k; for (i__ = lp1; i__ <= i__2; ++i__) { t = f01lzz_(&f, &x, &small, &big); f01lzw_(&t, &c__, &s, &sqteps, &rsqtps, &big); if (i__ > lp1) { work1[i__ - 1] = c__ * f + s * x; } f = c__ * sv[i__ - 1] + s * work1[i__]; work1[i__] = c__ * work1[i__] - s * sv[i__ - 1]; x = s * sv[i__]; svi = c__ * sv[i__]; if (! (*wantz)) { goto L100; } work2[i__] = c__; work3[i__] = s; L100: t = f01lzz_(&f, &x, &small, &big); f01lzw_(&t, &c__, &s, &sqteps, &rsqtps, &big); if (*wanty) { f01lzy_(ly, &c__, &s, &y[(i__ - 1) * y_dim1 + 1], &y[i__ * y_dim1 + 1]); } if (! (*wantb)) { goto L120; } t = b[i__]; b[i__] = c__ * t - s * b[i__ - 1]; b[i__ - 1] = c__ * b[i__ - 1] + s * t; L120: sv[i__ - 1] = c__ * f + s * x; f = c__ * work1[i__] + s * svi; sv[i__] = c__ * svi - s * work1[i__]; if (i__ == k) { goto L140; } x = s * work1[i__ + 1]; work1[i__ + 1] = c__ * work1[i__ + 1]; L140: ; } work1[k] = f; if (! (*wantz)) { goto L40; } i__2 = *ncz; for (j = 1; j <= i__2; ++j) { i__3 = k - l + 1; f02szz_(&i__3, &work2[l], &work3[l], &z__[l + j * z_dim1]); /* L160: */ } goto L40; /* COME TO NEXT PIECE IF SV(L-1) IS NEGLIGIBLE. FORCE A SPLIT. */ L180: lm1 = l; ++l; x = work1[l]; work1[l] = 0.; i__2 = k; for (i__ = l; i__ <= i__2; ++i__) { t = f01lzz_(&sv[i__], &x, &small, &big); f01lzw_(&t, &c__, &s, &sqteps, &rsqtps, &big); if (*wanty) { d__1 = -s; f01lzy_(ly, &c__, &d__1, &y[lm1 * y_dim1 + 1], &y[i__ * y_dim1 + 1]); } if (! (*wantb)) { goto L200; } t = b[i__]; b[i__] = c__ * t + s * b[lm1]; b[lm1] = c__ * b[lm1] - s * t; L200: sv[i__] = c__ * sv[i__] + s * x; if (i__ == k) { goto L220; } x = -s * work1[i__ + 1]; work1[i__ + 1] = c__ * work1[i__ + 1]; L220: ; } /* IF WE COME HERE WITH L=K THEN A SINGULAR VALUE HAS BEEN */ /* FOUND. */ L240: if (l < k) { goto L80; } --k; /* L260: */ } L280: *ifail = k - 1; work1[1] = (doublereal) iter; /* NOW MAKE SINGULAR VALUES NON-NEGATIVE. */ /* K WILL BE 1 UNLESS FAILURE HAS OCCURED. */ i__1 = *n; for (j = k; j <= i__1; ++j) { if (sv[j] >= 0.) { goto L320; } sv[j] = -sv[j]; if (*wantb) { b[j] = -b[j]; } if (! (*wanty)) { goto L320; } i__2 = *ly; for (i__ = 1; i__ <= i__2; ++i__) { y[i__ + j * y_dim1] = -y[i__ + j * y_dim1]; /* L300: */ } L320: ; } /* NOW SORT THE SINGULAR VALUES INTO DESCENDING ORDER. */ if (*wantz) { jj = 0; } i__1 = *n; for (j = k; j <= i__1; ++j) { s = 0.; l = j; i__2 = *n; for (i__ = j; i__ <= i__2; ++i__) { if (sv[i__] <= s) { goto L340; } s = sv[i__]; l = i__; L340: ; } if (s == 0.) { goto L420; } if (*wantz) { work2[j] = (doublereal) l; } if (l == j) { goto L400; } if (*wantz) { jj = j; } sv[l] = sv[j]; sv[j] = s; if (! (*wanty)) { goto L380; } i__2 = *ly; for (i__ = 1; i__ <= i__2; ++i__) { t = y[i__ + j * y_dim1]; y[i__ + j * y_dim1] = y[i__ + l * y_dim1]; y[i__ + l * y_dim1] = t; /* L360: */ } L380: if (! (*wantb)) { goto L400; } t = b[j]; b[j] = b[l]; b[l] = t; L400: ; } L420: if (! (*wantz)) { goto L480; } if (jj == 0) { goto L480; } i__1 = *ncz; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = jj; for (j = k; j <= i__2; ++j) { l = (integer) work2[j]; if (j == l) { goto L440; } t = z__[j + i__ * z_dim1]; z__[j + i__ * z_dim1] = z__[l + i__ * z_dim1]; z__[l + i__ * z_dim1] = t; L440: ; } /* L460: */ } L480: if (*ifail == 0) { return 0; } ++(*ifail); L500: *ifail = p01abf_(&ierr, ifail, "F02SZF", &c__0, p01rec, 6L, 1L); return 0; } /* f02szf_ */ /* Subroutine */ int f02szz_(n, c__, s, x) integer *n; doublereal *c__, *s, *x; { /* System generated locals */ integer i__1; /* Local variables */ static integer i__; static doublereal w; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (PLRT10) */ /* F02SZZ RETURNS THE N ELEMENT VECTOR */ /* Y = R(N-1,N)*R(N-2,N-1)*...*R(1,2)*X , */ /* WHERE X IS AN N ELEMENT VECTOR AND R(J-1,J) IS A PLANE */ /* ROTATION FOR THE (J-1,J)-PLANE. */ /* Y IS OVERWRITTEN ON X. */ /* THE N ELEMENT VECTORS C AND S MUST BE SUCH THAT THE */ /* NON-IDENTITY PART OF R(J-1,J) IS GIVEN BY */ /* R(J-1,J) = ( C(J) S(J) ) . */ /* ( -S(J) C(J) ) */ /* C(1) AND S(1) ARE NOT REFERENCED. */ /* N MUST BE AT LEAST 1. IF N=1 THEN AN IMMEDIATE RETURN TO */ /* THE CALLING PROGRAM IS MADE. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --x; --s; --c__; /* Function Body */ if (*n == 1) { return 0; } i__1 = *n; for (i__ = 2; i__ <= i__1; ++i__) { w = x[i__ - 1]; x[i__ - 1] = c__[i__] * w + s[i__] * x[i__]; x[i__] = c__[i__] * x[i__] - s[i__] * w; /* L20: */ } return 0; } /* f02szz_ */ /* Subroutine */ int f02waf_(m, n, a, nra, wantb, b, sv, work, lwork, ifail) integer *m, *n; doublereal *a; integer *nra; logical *wantb; doublereal *b, *sv, *work; integer *lwork, *ifail; { /* System generated locals */ integer a_dim1, a_offset; /* Local variables */ static integer ierr, npnp1; extern integer p01abf_(); extern /* Subroutine */ int f01qcf_(), f01qdf_(); static char p01rec[1*1]; extern /* Subroutine */ int f01lzf_(), f02way_(), f02szf_(); static integer np1; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 15 REVISED. IER-912 (APR 1991). */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (SVDGN1) */ /* Modified by Sven to replace calls to F01QAF and F02WAZ. */ /* 7-Feb-1991. */ /* F02WAF RETURNS PART OF THE SINGULAR VALUE DECOMPOSITION */ /* OF THE M*N (M.GE.N) MATRIX A GIVEN BY */ /* A = Q*(D)*(P**T) , */ /* (0) */ /* WHERE Q AND P ARE ORTHOGONAL MATRICES AND D IS AN N*N */ /* DIAGONAL MATRIX WITH NON-NEGATIVE DIAGONAL ELEMENTS, */ /* THESE BEING THE SINGULAR VALUES OF A. */ /* P**T AND THE DIAGONAL ELEMENTS OF D ARE RETURNED. */ /* IF WANTB IS .TRUE. THEN (Q**T)*B IS ALSO RETURNED. */ /* INPUT PARAMETERS. */ /* M - NUMBER OF ROWS OF A. M MUST BE AT LEAST N. */ /* N - NUMBER OF COLUMNS OF A. N MUST BE AT LEAST */ /* UNITY AND MUST NOT BE LARGER THAN THAN M. */ /* A - THE M*N MATRIX TO BE FACTORIZED. */ /* NRA - ROW DIMENSION OF A AS DECLARED IN THE CALLING PROGRAM. */ /* NRA MUST BE AT LEAST M. */ /* WANTB - MUST BE .TRUE. IF (Q**T)*B IS REQUIRED. */ /* IF WANTB IS .FALSE. THEN B IS NOT REFERENCED. */ /* B - AN M ELEMENT VECTOR. */ /* IFAIL - THE USUAL FAILURE PARAMETER. IF IN DOUBT SET */ /* IFAIL TO ZERO BEFORE CALLING F02WAF. */ /* OUTPUT PARAMETERS. */ /* A - THE TOP N*N PART OF A WILL CONTAIN THE N*N ORTHOGONAL */ /* MATRIX P**T. */ /* THE REMAINING (M-N)*N PART OF A IS USED FOR INTERNAL */ /* WORKSPACE. */ /* B - IF WANTB IS .TRUE. THEN B IS OVERWRITTEN BY */ /* THE M ELEMENT VECTOR (Q**T)*B. */ /* SV - N ELEMENT VECTOR CONTAINING THE SINGULAR */ /* VALUES OF A. THEY ARE ORDERED SO THAT */ /* SV(1).GE.SV(2).GE. ... .GE.SV(N).GE.0. */ /* IFAIL - ON NORMAL RETURN IFAIL WILL BE ZERO. */ /* IN THE UNLIKELY EVENT THAT THE QR-ALGORITHM */ /* FAILS TO FIND THE SINGULAR VALUES IN 50*N */ /* ITERATIONS THEN IFAIL IS SET TO 2. */ /* IF AN INPUT PARAMETER IS INCORRECTLY SUPPLIED */ /* THEN IFAIL IS SET TO UNITY. */ /* WORKSPACE PARAMETERS. */ /* WORK - A 3*N ELEMENT VECTOR. */ /* WORK(1) RETURNS THE TOTAL NUMBER OF ITERATIONS TAKEN */ /* BY THE QR-ALGORITHM. */ /* LWORK - THE LENGTH OF THE VECTOR WORK. LWORK MUST BE */ /* AT LEAST 3*N. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --b; --sv; a_dim1 = *nra; a_offset = a_dim1 + 1; a -= a_offset; --work; /* Function Body */ ierr = *ifail; if (ierr == 0) { *ifail = 1; } if (*nra < *m || *m < *n || *lwork < *n * 3 || *n < 1) { goto L20; } np1 = *n + 1; npnp1 = *n + np1; /* Call to F01QAF replaced by a call to F01QCF. */ /* CALL F01QAF(M,N,A,NRA,A,NRA,WORK,IFAIL) */ f01qcf_(m, n, &a[a_offset], nra, &work[1], ifail); /* Call to F02WAZ replaced by a call to F01QDF. */ /* IF (WANTB) CALL F02WAZ(M,N,A,NRA,WORK,B,B) */ if (*wantb) { f01qdf_("Transpose", "Separate", m, n, &a[a_offset], nra, &work[1], & c__1, &b[1], m, &work[*n + 1], ifail, 9L, 8L); } f01lzf_(n, &a[a_offset], nra, &a[a_offset], nra, wantb, &b[1], &c_false, & c_false, &work[1], &c__1, &c__1, &c_false, &work[1], &c__1, &c__1, &sv[1], &work[1], &work[1], &work[1], ifail); f02way_(n, &a[a_offset], nra, &a[a_offset], nra); *ifail = 1; f02szf_(n, &sv[1], &work[1], &sv[1], wantb, &b[1], &c_false, &work[1], & c__1, &c__1, &c_true, &a[a_offset], nra, n, &work[1], &work[np1], &work[npnp1], ifail); if (*ifail == 0) { return 0; } *ifail = 2; L20: *ifail = p01abf_(&ierr, ifail, "F02WAF", &c__0, p01rec, 6L, 1L); return 0; } /* f02waf_ */ /* Subroutine */ int f02way_(n, c__, nrc, pt, nrpt) integer *n; doublereal *c__; integer *nrc; doublereal *pt; integer *nrpt; { /* System generated locals */ integer c_dim1, c_offset, pt_dim1, pt_offset, i__1, i__2, i__3; doublereal d__1; /* Builtin functions */ double sqrt(); /* Local variables */ static integer i__, j, k; static doublereal t; extern doublereal x02ajf_(), x02amf_(); extern /* Subroutine */ int f01lzw_(), f01lzy_(); static doublereal cs; static integer kk; static doublereal sn; static integer km1, kp1; static doublereal sqteps, rsqtps, big; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 12 REVISED. IER-519 (AUG 1986). */ /* MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988). */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (BIGVPT) */ /* F02WAY RETURNS THE N*N ORTHOGONAL MATRIX P**T FOR THE */ /* FACTORIZATION OF ROUTINE F01LZF. */ /* DETAILS OF P MUST BE SUPPLIED IN THE N*N MATRIX C AS */ /* RETURNED FROM ROUTINE F01LZF. */ /* NRC AND NRPT MUST BE THE ROW DIMENSIONS OF C AND PT */ /* RESPECTIVELY AS DECLARED IN THE CALLING PROGRAM AND MUST */ /* EACH BE AT LEAST N. */ /* THE ROUTINE MAY BE CALLED WITH PT=C. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ c_dim1 = *nrc; c_offset = c_dim1 + 1; c__ -= c_offset; pt_dim1 = *nrpt; pt_offset = pt_dim1 + 1; pt -= pt_offset; /* Function Body */ big = 1. / x02amf_(); sqteps = sqrt(x02ajf_()); rsqtps = 1. / sqteps; i__1 = *n; for (j = 3; j <= i__1; ++j) { i__2 = j - 2; for (i__ = 1; i__ <= i__2; ++i__) { pt[i__ + j * pt_dim1] = c__[i__ + j * c_dim1]; /* L5: */ } /* L15: */ } pt[*n + *n * pt_dim1] = 1.; if (*n == 1) { return 0; } pt[*n - 1 + *n * pt_dim1] = 0.; pt[*n + (*n - 1) * pt_dim1] = 0.; pt[*n - 1 + (*n - 1) * pt_dim1] = 1.; if (*n == 2) { return 0; } k = *n; i__1 = *n; for (kk = 3; kk <= i__1; ++kk) { kp1 = k; --k; km1 = k - 1; pt[km1 + k * pt_dim1] = 0.; i__2 = *n; for (j = kp1; j <= i__2; ++j) { t = pt[km1 + j * pt_dim1]; pt[km1 + j * pt_dim1] = 0.; if (t == 0.) { goto L20; } d__1 = -t; f01lzw_(&d__1, &cs, &sn, &sqteps, &rsqtps, &big); i__3 = *n - km1; f01lzy_(&i__3, &cs, &sn, &pt[k + (j - 1) * pt_dim1], &pt[k + j * pt_dim1]); L20: ; } pt[km1 + km1 * pt_dim1] = 1.; i__2 = *n; for (i__ = k; i__ <= i__2; ++i__) { pt[i__ + km1 * pt_dim1] = 0.; /* L40: */ } /* L60: */ } return 0; } /* f02way_ */ integer f02wdy_(n, sv, tol) integer *n; doublereal *sv, *tol; { /* System generated locals */ integer ret_val, i__1; /* Local variables */ static integer i__; static doublereal delta; extern doublereal x02ajf_(); static integer ir; static doublereal tl; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988). */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (IRKSVD) */ /* F02WDY RETURNS THE RANK OF AN M*K MATRIX A FOLLOWING A */ /* SINGULAR VALUE DECOMPOSITION OF A. */ /* THE N=MIN(M,K) SINGULAR VALUES OF A MUST BE IN */ /* DESCENDING ORDER IN THE N ELEMENT VECTOR SV. THEN F02WDY */ /* RETURNS THE LARGEST INTEGER SUCH THAT */ /* SV(F02WDY) .GT. TOL*SV(1) . */ /* IF SV(1)=0 THEN F02WDY IS RETURNED AS ZERO. */ /* IF TOL.LT.EPS OR TOL.GE.1 THEN THE VALUE EPS IS USED IN */ /* PLACE OF TOL, WHERE EPS IS THE SMALLEST REAL FOR WHICH */ /* 1.0+EPS.GT.1.0 ON THE MACHINE. FOR MOST PROBLEMS THIS IS */ /* UNREASONABLY SMALL AND TOL SHOULD BE CHOSEN TO */ /* APPROXIMATE THE RELATIVE ERRORS IN THE ELEMENTS OF A. */ /* IF INSTEAD SINGULAR VALUES BELOW SOME VALUE DELTA ARE TO */ /* BE REGARDED AS ZERO THEN SUPPLY TOL AS DELTA/SV(1). */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --sv; /* Function Body */ tl = *tol; delta = x02ajf_(); if (tl < delta || tl >= 1.) { tl = delta; } ir = 0; delta = tl * sv[1]; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (sv[i__] <= delta) { goto L40; } ir = i__; /* L20: */ } L40: ret_val = ir; return ret_val; } /* f02wdy_ */ /* Subroutine */ int f04aef_(a, ia, b, ib, n, m, c__, ic, wkspce, aa, iaa, bb, ibb, ifail) doublereal *a; integer *ia; doublereal *b; integer *ib, *n, *m; doublereal *c__; integer *ic; doublereal *wkspce, *aa; integer *iaa; doublereal *bb; integer *ibb, *ifail; { /* Format strings */ static char fmt_99999[] = "(1x,\002** On entry, N.lt.0: N =\002,i16)"; static char fmt_99998[] = "(1x,\002** On entry, M.lt.0: M =\002,i16)"; static char fmt_99997[] = "(1x,\002** On entry, IA.lt.max(1,N): IA =\002\ ,i16,\002, N =\002,i16)"; static char fmt_99996[] = "(1x,\002** On entry, IB.lt.max(1,N): IB =\002\ ,i16,\002, N =\002,i16)"; static char fmt_99995[] = "(1x,\002** On entry, IC.lt.max(1,N): IC =\002\ ,i16,\002, N =\002,i16)"; static char fmt_99994[] = "(1x,\002** On entry, IAA.lt.max(1,N): IAA \ =\002,i16,\002, N =\002,i16)"; static char fmt_99993[] = "(1x,\002** On entry, IBB.lt.max(1,N): IBB \ =\002,i16,\002, N =\002,i16)"; static char fmt_99991[] = "(1x,\002** Matrix A is too ill-conditioned\ .\002)"; static char fmt_99992[] = "(1x,\002** Matrix A is approximately singul\ ar.\002)"; /* System generated locals */ integer a_dim1, a_offset, aa_dim1, aa_offset, b_dim1, b_offset, bb_dim1, bb_offset, c_dim1, c_offset, i__1, i__2; doublereal d__1; /* Builtin functions */ integer s_wsfi(), do_fio(), e_wsfi(); /* Local variables */ static integer nrec, info, ierr; extern /* Subroutine */ int f04ahf_(); static integer i__, j; extern integer p01abf_(); extern /* Subroutine */ int f07adg_(); static char p01rec[80*1]; extern doublereal x02ajf_(); static integer ifail1; /* Fortran I/O blocks */ static icilist io___198 = { 0, p01rec, 0, fmt_99999, 80, 1 }; static icilist io___199 = { 0, p01rec, 0, fmt_99998, 80, 1 }; static icilist io___200 = { 0, p01rec, 0, fmt_99997, 80, 1 }; static icilist io___201 = { 0, p01rec, 0, fmt_99996, 80, 1 }; static icilist io___202 = { 0, p01rec, 0, fmt_99995, 80, 1 }; static icilist io___203 = { 0, p01rec, 0, fmt_99994, 80, 1 }; static icilist io___204 = { 0, p01rec, 0, fmt_99993, 80, 1 }; static icilist io___209 = { 0, p01rec, 0, fmt_99991, 80, 1 }; static icilist io___210 = { 0, p01rec, 0, fmt_99992, 80, 1 }; /* MARK 15 RE-ISSUE. NAG COPYRIGHT 1991. */ /* Accurate solution of a set of real linear equations */ /* with multiple right hand sides. */ /* 1st August 1971 */ /* Rewritten to call F07ADG, a modified version of LAPACK routine */ /* SGETRF/F07ADF; new IFAIL exit inserted for illegal input */ /* parameters; error messages inserted. February 1991. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ a_dim1 = *ia; a_offset = a_dim1 + 1; a -= a_offset; b_dim1 = *ib; b_offset = b_dim1 + 1; b -= b_offset; c_dim1 = *ic; c_offset = c_dim1 + 1; c__ -= c_offset; --wkspce; aa_dim1 = *iaa; aa_offset = aa_dim1 + 1; aa -= aa_offset; bb_dim1 = *ibb; bb_offset = bb_dim1 + 1; bb -= bb_offset; /* Function Body */ ierr = 0; nrec = 0; if (*n < 0) { ierr = 3; nrec = 1; s_wsfi(&io___198); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfi(); } else if (*m < 0) { ierr = 3; nrec = 1; s_wsfi(&io___199); do_fio(&c__1, (char *)&(*m), (ftnlen)sizeof(integer)); e_wsfi(); } else if (*ia < max(1,*n)) { ierr = 3; nrec = 1; s_wsfi(&io___200); do_fio(&c__1, (char *)&(*ia), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfi(); } else if (*ib < max(1,*n)) { ierr = 3; nrec = 1; s_wsfi(&io___201); do_fio(&c__1, (char *)&(*ib), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfi(); } else if (*ic < max(1,*n)) { ierr = 3; nrec = 1; s_wsfi(&io___202); do_fio(&c__1, (char *)&(*ic), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfi(); } else if (*iaa < max(1,*n)) { ierr = 3; nrec = 1; s_wsfi(&io___203); do_fio(&c__1, (char *)&(*iaa), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfi(); } else if (*ibb < max(1,*n)) { ierr = 3; nrec = 1; s_wsfi(&io___204); do_fio(&c__1, (char *)&(*ibb), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfi(); } else { /* Copy A to working array AA */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *n; for (j = 1; j <= i__2; ++j) { aa[i__ + j * aa_dim1] = a[i__ + j * a_dim1]; /* L20: */ } /* L40: */ } f07adg_(n, n, &aa[aa_offset], iaa, &wkspce[1], &info); if (info == 0) { if (*n > 0 && *m > 0) { ifail1 = 1; d__1 = x02ajf_(); f04ahf_(n, m, &a[a_offset], ia, &aa[aa_offset], iaa, &wkspce[ 1], &b[b_offset], ib, &d__1, &c__[c_offset], ic, &bb[ bb_offset], ibb, &i__, &ifail1); if (ifail1 != 0) { ierr = 2; nrec = 1; s_wsfi(&io___209); e_wsfi(); } } } else { ierr = 1; nrec = 1; s_wsfi(&io___210); e_wsfi(); } } *ifail = p01abf_(ifail, &ierr, "F04AEF", &nrec, p01rec, 6L, 80L); return 0; } /* f04aef_ */ /* Subroutine */ int f04ahf_(n, ir, a, ia, aa, iaa, p, b, ib, eps, x, ix, bb, ibb, l, ifail) integer *n, *ir; doublereal *a; integer *ia; doublereal *aa; integer *iaa; doublereal *p, *b; integer *ib; doublereal *eps, *x; integer *ix; doublereal *bb; integer *ibb, *l, *ifail; { /* System generated locals */ integer a_dim1, a_offset, aa_dim1, aa_offset, b_dim1, b_offset, bb_dim1, bb_offset, x_dim1, x_offset, i__1, i__2, i__3, i__4; doublereal d__1, d__2; /* Local variables */ static doublereal xmax; static integer i__, j; extern integer p01abf_(); extern /* Subroutine */ int f04ajf_(), x03aaf_(); static doublereal bbmax; static char p01rec[1*1]; static integer isave; static doublereal d0, d1, d2; static integer ifail1; static doublereal d11; static integer id2; /* MARK 2 RELEASE. NAG COPYRIGHT 1972 */ /* MARK 3 REVISED. */ /* MARK 4 REVISED. */ /* MARK 4.5 REVISED */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* UNSYMACCSOLVE */ /* SOLVES AX=B WHERE A IS AN N*N UNSYMMETRIC MATRIX AND B IS AN */ /* N*IR MATRIX OF RIGHT HAND SIDES, USING THE SUBROUTINE F04AJF. */ /* THE SUBROUTINE MUST BY PRECEDED BY F03AFF IN WHICH L AND U */ /* ARE PRODUCED IN AA(I,J) AND THE INTERCHANGES IN P(I). THE */ /* RESIDUALS BB=B-AX ARE CALCULATED AND AD=BB IS SOLVED, OVER- */ /* WRITING D ON BB. THE REFINEMENT IS REPEATED, AS LONG AS THE */ /* MAXIMUM CORRECTION AT ANY STAGE IS LESS THAN HALF THAT AT THE */ /* PREVIOUS STAGE, UNTIL THE MAXIMUM CORRECTION IS LESS THAN 2 */ /* EPS TIMES THE MAXIMUM X. SETS IFAIL = 1 IF THE SOLUTION FAILS */ /* TO IMPROVE, ELSE IFAIL = 0. L IS THE NUMBER OF ITERATIONS. */ /* ADDITIONAL PRECISION INNERPRODUCTS ARE ABSOLUTELY NECESSARY. */ /* 1ST DECEMBER 1971 */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --p; a_dim1 = *ia; a_offset = a_dim1 + 1; a -= a_offset; aa_dim1 = *iaa; aa_offset = aa_dim1 + 1; aa -= aa_offset; b_dim1 = *ib; b_offset = b_dim1 + 1; b -= b_offset; x_dim1 = *ix; x_offset = x_dim1 + 1; x -= x_offset; bb_dim1 = *ibb; bb_offset = bb_dim1 + 1; bb -= bb_offset; /* Function Body */ isave = *ifail; ifail1 = 0; i__1 = *ir; for (j = 1; j <= i__1; ++j) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { x[i__ + j * x_dim1] = 0.; bb[i__ + j * bb_dim1] = b[i__ + j * b_dim1]; /* L20: */ } /* L40: */ } *l = 0; d0 = 0.; L60: f04ajf_(n, ir, &aa[aa_offset], iaa, &p[1], &bb[bb_offset], ibb); ++(*l); id2 = 0; d1 = 0.; i__1 = *ir; for (j = 1; j <= i__1; ++j) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { x[i__ + j * x_dim1] += bb[i__ + j * bb_dim1]; /* L80: */ } /* L100: */ } i__1 = *ir; for (j = 1; j <= i__1; ++j) { xmax = 0.; bbmax = 0.; i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { if ((d__1 = x[i__ + j * x_dim1], abs(d__1)) > xmax) { xmax = (d__2 = x[i__ + j * x_dim1], abs(d__2)); } if ((d__1 = bb[i__ + j * bb_dim1], abs(d__1)) > bbmax) { bbmax = (d__2 = bb[i__ + j * bb_dim1], abs(d__2)); } i__3 = *n * *ia - i__ + 1; i__4 = (*ir - j + 1) * *ix; d__1 = -b[i__ + j * b_dim1]; x03aaf_(&a[i__ + a_dim1], &i__3, &x[j * x_dim1 + 1], &i__4, n, ia, &c__1, &d__1, &c_b114, &d11, &d2, &c_true, &ifail1); bb[i__ + j * bb_dim1] = -d11; /* L120: */ } if (bbmax > d1 * xmax) { d1 = bbmax / xmax; } if (bbmax > *eps * 2. * xmax) { id2 = 1; } /* L140: */ } if (d1 > d0 * .5 && *l != 1) { goto L160; } d0 = d1; if (id2 == 1) { goto L60; } *ifail = 0; return 0; L160: *ifail = p01abf_(&isave, &c__1, "F04AHF", &c__0, p01rec, 6L, 1L); return 0; } /* f04ahf_ */ /* Subroutine */ int f04ajf_(n, ir, a, ia, p, b, ib) integer *n, *ir; doublereal *a; integer *ia; doublereal *p, *b; integer *ib; { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2; /* Local variables */ static integer i__, k; static doublereal x; static integer i1; extern /* Subroutine */ int dtrsv_(); /* MARK 2 RELEASE. NAG COPYRIGHT 1972 */ /* MARK 4 REVISED. */ /* MARK 4.5 REVISED */ /* MARK 11 REVISED. VECTORISATION (JAN 1984). */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 12 REVISED. EXTENDED BLAS (JUNE 1986) */ /* UNSYMSOL */ /* SOLVES AX=B, WHERE A IS AN UNSYMMETRIC MATRIX AND B IS AN */ /* N*IR */ /* MATRIX OF IR RIGHT-HAND SIDES. THE SUBROUTINE F04AJF MUST BY */ /* PRECEDED BY F03AFF IN WHICH L AND U ARE PRODUCED IN A(I,J), */ /* FROM A, AND THE RECORD OF THE INTERCHANGES IS PRODUCED IN */ /* P(I). AX=B IS SOLVED IN THREE STEPS, INTERCHANGE THE */ /* ELEMENTS OF B, LY=B AND UX=Y. THE MATRICES Y AND THEN X ARE */ /* OVERWRITTEN ON B. */ /* 1ST AUGUST 1971 */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Executable Statements .. */ /* INTERCHANGING ELEMENTS OF B */ /* Parameter adjustments */ --p; a_dim1 = *ia; a_offset = a_dim1 + 1; a -= a_offset; b_dim1 = *ib; b_offset = b_dim1 + 1; b -= b_offset; /* Function Body */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i1 = (integer) (p[i__] + .5); if (i1 == i__) { goto L40; } i__2 = *ir; for (k = 1; k <= i__2; ++k) { x = b[i__ + k * b_dim1]; b[i__ + k * b_dim1] = b[i1 + k * b_dim1]; b[i1 + k * b_dim1] = x; /* L20: */ } L40: ; } i__1 = *ir; for (k = 1; k <= i__1; ++k) { /* SOLUTION OF LY= B */ dtrsv_("L", "N", "N", n, &a[a_offset], ia, &b[k * b_dim1 + 1], &c__1, 1L, 1L, 1L); /* SOLUTION OF UX= Y */ dtrsv_("U", "N", "U", n, &a[a_offset], ia, &b[k * b_dim1 + 1], &c__1, 1L, 1L, 1L); /* L60: */ } return 0; } /* f04ajf_ */ /* Subroutine */ int f04jaf_(m, n, a, nra, b, tol, sigma, irank, work, lwork, ifail) integer *m, *n; doublereal *a; integer *nra; doublereal *b, *tol, *sigma; integer *irank; doublereal *work; integer *lwork, *ifail; { /* System generated locals */ integer a_dim1, a_offset; /* Local variables */ static integer ierr; extern integer p01abf_(); extern /* Subroutine */ int f02waf_(); static char p01rec[1*1]; extern /* Subroutine */ int f04jaz_(); extern integer f02wdy_(); static integer np1, np2, nnn; /* MARK 8 RELEASE. NAG COPYRIGHT 1979. */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (SVDLS1) */ /* F04JAF RETURNS THE N ELEMENT VECTOR X, OF MINIMAL */ /* LENGTH, THAT MINIMIZES THE EUCLIDEAN LENGTH OF THE M */ /* ELEMENT VECTOR R GIVEN BY */ /* R = B-A*X , */ /* WHERE A IS AN M*N (M.GE.N) MATRIX AND B IS AN M ELEMENT */ /* VECTOR. X IS OVERWRITTEN ON B. */ /* THE SOLUTION IS OBTAINED VIA A SINGULAR VALUE */ /* DECOMPOSITION (SVD) OF THE MATRIX A GIVEN BY */ /* A = Q*(D)*(P**T) , */ /* (0) */ /* WHERE Q AND P ARE ORTHOGONAL AND D IS A DIAGONAL MATRIX WITH */ /* NON-NEGATIVE DIAGONAL ELEMENTS, THESE BEING THE SINGULAR */ /* VALUES OF A. */ /* INPUT PARAMETERS. */ /* M - NUMBER OF ROWS OF A. M MUST BE AT LEAST N. */ /* N - NUMBER OF COLUMNS OF A. N MUST BE AT LEAST UNITY. */ /* A - AN M*N REAL MATRIX. */ /* NRA - ROW DIMENSION OF A AS DECLARED IN THE CALLING PROGRAM. */ /* NRA MUST BE AT LEAST M. */ /* B - AN M ELEMENT REAL VECTOR. */ /* TOL - A RELATIVE TOLERANCE USED TO DETERMINE THE RANK OF A. */ /* TOL SHOULD BE CHOSEN AS APPROXIMATELY THE */ /* LARGEST RELATIVE ERROR IN THE ELEMENTS OF A. */ /* FOR EXAMPLE IF THE ELEMENTS OF A ARE CORRECT */ /* TO ABOUT 4 SIGNIFICANT FIGURES THEN TOL */ /* SHOULD BE CHOSEN AS ABOUT 5.0*10.0**(-4). */ /* IFAIL - THE USUAL FAILURE PARAMETER. IF IN DOUBT SET */ /* IFAIL TO ZERO BEFORE CALLING THIS ROUTINE. */ /* OUTPUT PARAMETERS. */ /* A - THE TOP N*N PART OF A WILL CONTAIN THE */ /* ORTHOGONAL MATRIX P**T OF THE SVD. */ /* THE REMAINDER OF A IS USED FOR INTERNAL WORKSPACE. */ /* B - THE FIRST N ELEMENTS OF B WILL CONTAIN THE */ /* MINIMAL LEAST SQUARES SOLUTION VECTOR X. */ /* SIGMA - IF M IS GREATER THAN IRANK THEN SIGMA WILL CONTAIN THE */ /* STANDARD ERROR GIVEN BY */ /* SIGMA=L(R)/SQRT(M-IRANK), WHERE L(R) DENOTES */ /* THE EUCLIDEAN LENGTH OF THE RESIDUAL VECTOR */ /* R. IF M=IRANK THEN SIGMA IS RETURNED AS ZERO. */ /* IRANK - THE RANK OF THE MATRIX A. */ /* IFAIL - ON NORMAL RETURN IFAIL WILL BE ZERO. */ /* IN THE UNLIKELY EVENT THAT THE QR-ALGORITHM */ /* FAILS TO FIND THE SINGULAR VALUES IN 50*N */ /* ITERATIONS THEN IFAIL IS SET TO 2. */ /* IF AN INPUT PARAMETER IS INCORRECTLY SUPPLIED */ /* THEN IFAIL IS SET TO UNITY. */ /* WORKSPACE PARAMETERS. */ /* WORK - A 4*N ELEMENT VECTOR. */ /* ON RETURN THE FIRST N ELEMENTS OF WORK WILL */ /* CONTAIN THE SINGULAR VALUES OF A ARRANGED IN */ /* DESCENDING ORDER. */ /* WORK(N+1) WILL CONTAIN THE TOTAL NUMBER OF ITERATIONS */ /* TAKEN BY THE QR-ALGORITHM. */ /* LWORK - THE LENGTH OF THE VECTOR WORK. LWORK MUST BE */ /* AT LEAST 4*N. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --b; a_dim1 = *nra; a_offset = a_dim1 + 1; a -= a_offset; --work; /* Function Body */ ierr = *ifail; if (ierr == 0) { *ifail = 1; } if (*nra < *m || *m < *n || *n < 1 || *lwork < *n << 2) { goto L20; } np1 = *n + 1; np2 = np1 + 1; nnn = *n * 3; f02waf_(m, n, &a[a_offset], nra, &c_true, &b[1], &work[1], &work[np1], & nnn, ifail); if (*ifail != 0) { goto L20; } *irank = f02wdy_(n, &work[1], tol); f04jaz_(m, n, irank, &work[1], n, &b[1], &a[a_offset], nra, &b[1], sigma, &work[np2]); return 0; L20: *ifail = p01abf_(&ierr, ifail, "F04JAF", &c__0, p01rec, 6L, 1L); return 0; } /* f04jaf_ */ /* Subroutine */ int f04jay_(n, irank, sv, lsv, b, pt, nrpt, x, work) integer *n, *irank; doublereal *sv; integer *lsv; doublereal *b, *pt; integer *nrpt; doublereal *x, *work; { /* System generated locals */ integer pt_dim1, pt_offset, i__1; /* Local variables */ static integer i__; extern /* Subroutine */ int dgemv_(); /* MARK 13 RE-ISSUE. NAG COPYRIGHT 1988. */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (SVDLSQ) */ /* F04JAY RETURNS THE N ELEMENT VECTOR X GIVEN BY */ /* X = P*(D**(-1))*B , */ /* WHERE D IS AN IRANK*IRANK NON-SINGULAR DIAGONAL MATRIX, */ /* P CONTAINS THE FIRST IRANK COLUMNS OF AN N*N ORTHOGONAL */ /* MATRIX AND B IS AN IRANK ELEMENT VECTOR. */ /* THE ROUTINE MAY BE CALLED WITH IRANK=0 IN WHICH CASE X */ /* IS RETURNED AS THE ZERO VECTOR. */ /* INPUT PARAMETERS. */ /* N - NUMBER OF ROWS OF P. N MUST BE AT LEAST UNITY. */ /* IRANK - ORDER OF THE MATRIX D. */ /* IF IRANK=0 THEN SV, B, PT AND WORK ARE NOT REFERENCED. */ /* SV - AN IRANK ELEMENT VECTOR CONTAINING THE */ /* DIAGONAL ELEMENTS OF D. SV MUST BE SUCH THAT */ /* NO ELEMENT OF (D**(-1)*B WILL OVERFLOW. */ /* LSV - LSV MUST BE AT LEAST MAX(1,IRANK). */ /* B - AN IRANK ELEMENT VECTOR. */ /* PT - AN IRANK*N ELEMENT MATRIX CONTAINING THE MATRIX P**T. */ /* NRPT - ROW DIMENSION OF PT AS DECLARED IN THE */ /* CALLING PROGRAM. NRPT MUST BE AT LEAST LSV. */ /* OUTPUT PARAMETER. */ /* X - N ELEMENT VECTOR CONTAINING P*(D**(-1))*B. */ /* IF IRANK=0 THEN X RETURNS THE ZERO VECTOR. */ /* THE ROUTINE MAY BE CALLED WITH X=B OR WITH X=SV. */ /* WORKSPACE PARAMETER. */ /* WORK - AN LSV ELEMENT VECTOR. */ /* IF THE ROUTINE IS NOT CALLED WITH X=B THEN IT MAY BE */ /* CALLED WITH WORK=B. SIMILARLY IF THE ROUTINE */ /* IS NOT CALLED WITH X=SV THEN IT MAY BE CALLED */ /* WITH WORK=SV. */ /* Modified to call BLAS. */ /* Jeremy Du Croz, NAG Central Office, October 1987. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --x; --work; --b; --sv; pt_dim1 = *nrpt; pt_offset = pt_dim1 + 1; pt -= pt_offset; /* Function Body */ if (*irank == 0) { goto L40; } i__1 = *irank; for (i__ = 1; i__ <= i__1; ++i__) { work[i__] = b[i__] / sv[i__]; /* L20: */ } dgemv_("Transpose", irank, n, &c_b15, &pt[pt_offset], nrpt, &work[1], & c__1, &c_b114, &x[1], &c__1, 9L); return 0; L40: i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { x[i__] = 0.; /* L60: */ } return 0; } /* f04jay_ */ /* Subroutine */ int f04jaz_(m, n, irank, sv, lsv, b, pt, nrpt, x, sigma, work) integer *m, *n, *irank; doublereal *sv; integer *lsv; doublereal *b, *pt; integer *nrpt; doublereal *x, *sigma, *work; { /* System generated locals */ integer pt_dim1, pt_offset; /* Builtin functions */ double sqrt(); /* Local variables */ static integer mmir; extern doublereal f06ejf_(); extern /* Subroutine */ int f04jay_(); static integer irp1; /* MARK 13 RE-ISSUE. NAG COPYRIGHT 1988. */ /* WRITTEN BY S. HAMMARLING, MIDDLESEX POLYTECHNIC (SVDLS0) */ /* F04JAZ RETURNS THE N ELEMENT VECTOR X, OF MINIMAL */ /* LENGTH, THAT MINIMIZES THE EUCLIDEAN LENGTH OF THE M */ /* ELEMENT VECTOR R GIVEN BY */ /* R = B-A*X , */ /* WHERE B IS AN M ELEMENT VECTOR AND A IS AN M*N MATRIX, */ /* FOLLOWING A SINGULAR VALUE DECOMPOSITION (SVD) OF A */ /* GIVEN BY */ /* A = Q*D*(P**T) , */ /* WHERE D IS A RECTANGULAR DIAGONAL MATRIX WHOSE DIAGONAL */ /* ELEMENTS CONTAIN THE SINGULAR VALUES OF A IN DESCENDING */ /* ORDER. */ /* INPUT PARAMETERS. */ /* M - NUMBER OF ROWS OF A. M MUST BE AT LEAST UNITY. */ /* N - NUMBER OF COLUMNS OF A. N MUST BE AT LEAST UNITY. */ /* IRANK - THE RANK OF THE MATRIX A. IRANK MUST BE SUCH THAT THE */ /* ELEMENTS SV(I), I=1,2,...,IRANK ARE NON-NEGLIGIBLE. */ /* IRANK MUST BE AT LEAST ZERO AND MUST NOT BE */ /* LARGER THAN MIN(M,N). */ /* ROUTINE F02WDY CAN BE USED TO DETERMINE RANK FOLLOWING */ /* AN SVD. */ /* SV - AN LSV ELEMENT VECTOR CONTAINING THE POSITIVE */ /* NON-NEGLIGIBLE SINGULAR VALUES OF A. */ /* LSV - LENGTH OF THE VECTOR SV. */ /* LSV MUST BE AT LEAST MAX(1,IRANK). */ /* B - MUST CONTAIN THE M ELEMENT VECTOR (Q**T)*B, WHERE Q IS */ /* THE LEFT-HAND ORTHOGONAL MATRIX OF THE SVD. */ /* PT - THE IRANK*N PART OF PT MUST CONTAIN THE FIRST */ /* IRANK ROWS OF THE RIGHT-HAND ORTHOGONAL */ /* MATRIX P**T OF THE SVD. */ /* NRPT - ROW DIMENSION OF PT AS DECLARED IN THE CALLING PROGRAM */ /* NRPT MUST BE AT LEAST LSV. */ /* OUTPUT PARAMETERS. */ /* X - THE N ELEMENT SOLUTION VECTOR. */ /* THE ROUTINE MAY BE CALLED WITH X=B OR WITH X=SV. */ /* SIGMA - IF M IS GREATER THAN IRANK THEN SIGMA WILL CONTAIN THE */ /* STANDARD ERROR GIVEN BY */ /* SIGMA=L(R)/SQRT(M-IRANK), WHERE L(R) DENOTES */ /* THE EUCLIDEAN LENGTH OF THE RESIDUAL VECTOR */ /* R. IF M=IRANK THEN SIGMA IS RETURNED AS ZERO. */ /* WORKSPACE PARAMETER. */ /* WORK - AN LSV ELEMENT VECTOR. */ /* IF THE ROUTINE IS NOT CALLED WITH X=B THEN IT MAY BE */ /* CALLED WITH WORK=B. SIMILARLY IF THE ROUTINE */ /* IS NOT CALLED WITH X=SV THEN IT MAY BE CALLED */ /* WITH WORK=SV. */ /* Modified to call BLAS. */ /* Jeremy Du Croz, NAG Central Office, October 1987. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --b; --x; --work; --sv; pt_dim1 = *nrpt; pt_offset = pt_dim1 + 1; pt -= pt_offset; /* Function Body */ *sigma = 0.; if (*irank == *m) { goto L20; } irp1 = *irank + 1; mmir = *m - *irank; *sigma = f06ejf_(&mmir, &b[irp1], &c__1) / sqrt((doublereal) mmir); L20: f04jay_(n, irank, &sv[1], lsv, &b[1], &pt[pt_offset], nrpt, &x[1], &work[ 1]); return 0; } /* f04jaz_ */ /* Subroutine */ int f06aaz_(srname, info, srname_len) char *srname; integer *info; ftnlen srname_len; { /* Format strings */ static char fmt_99999[] = "(\002 ** On entry to \002,a13,\002 parameter \ number \002,i2,\002 had an illegal value\002)"; /* Builtin functions */ integer s_wsfi(), do_fio(), e_wsfi(), s_cmp(); /* Subroutine */ int s_copy(); /* Local variables */ static integer ierr; extern integer p01acf_(); static integer ifail; static char varbnm[4], rec[80*1]; /* Fortran I/O blocks */ static icilist io___236 = { 0, rec, 0, fmt_99999, 80, 1 }; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* MARK 15 REVISED. IER-915 (APR 1991). */ /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* F06AAZ is an error handler for the Level 2 BLAS routines. */ /* It is called by the Level 2 BLAS routines if an input parameter is */ /* invalid. */ /* Parameters */ /* ========== */ /* SRNAME - CHARACTER*13. */ /* On entry, SRNAME specifies the name of the routine which */ /* called F06AAZ. */ /* INFO - INTEGER. */ /* On entry, INFO specifies the position of the invalid */ /* parameter in the parameter-list of the calling routine. */ /* Auxiliary routine for Level 2 Blas. */ /* Written on 20-July-1986. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. */ /* .. Executable Statements .. */ s_wsfi(&io___236); do_fio(&c__1, srname, 13L); do_fio(&c__1, (char *)&(*info), (ftnlen)sizeof(integer)); e_wsfi(); if (s_cmp(srname, "F06", 3L, 3L) == 0) { ierr = -1; s_copy(varbnm, " ", 4L, 4L); } else { ierr = -(*info); s_copy(varbnm, "INFO", 4L, 4L); } ifail = 0; ifail = p01acf_(&ifail, &ierr, srname, varbnm, &c__1, rec, 6L, 4L, 80L); return 0; /* End of F06AAZ. */ } /* f06aaz_ */ doublereal f06bmf_(scale, ssq) doublereal *scale, *ssq; { /* Initialized data */ static logical first = TRUE_; /* System generated locals */ doublereal ret_val; /* Builtin functions */ double sqrt(); /* Local variables */ static doublereal norm; extern doublereal x02amf_(); static doublereal flmin, flmax, sqt; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* .. Scalar Arguments .. */ /* .. */ /* F06BMF returns the value norm given by */ /* norm = ( scale*sqrt( ssq ), scale*sqrt( ssq ) .lt. flmax */ /* ( */ /* ( flmax, scale*sqrt( ssq ) .ge. flmax */ /* via the function name. */ /* Nag Fortran 77 O( 1 ) basic linear algebra routine. */ /* -- Written on 22-October-1982. */ /* Sven Hammarling, Nag Central Office. */ /* .. Local Scalars .. */ /* .. External Functions .. */ /* .. Intrinsic Functions .. */ /* .. Save statement .. */ /* .. Data statements .. */ /* .. */ /* .. Executable Statements .. */ if (first) { first = FALSE_; flmin = x02amf_(); flmax = 1 / flmin; } sqt = sqrt(*ssq); if (*scale < flmax / sqt) { norm = *scale * sqt; } else { norm = flmax; } ret_val = norm; return ret_val; /* End of F06BMF. ( SNORM ) */ } /* f06bmf_ */ /* Subroutine */ int f06edf_0_(n__, n, alpha, x, incx) int n__; integer *n; doublereal *alpha, *x; integer *incx; { /* System generated locals */ integer i__1, i__2; /* Local variables */ static integer ix; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* .. Entry Points .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* F06EDF performs the operation */ /* x := alpha*x */ /* Nag Fortran 77 version of the Blas routine DSCAL. */ /* Nag Fortran 77 O( n ) basic linear algebra routine. */ /* -- Written on 26-November-1982. */ /* Sven Hammarling, Nag Central Office. */ /* .. Parameters .. */ /* .. Local Scalars .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --x; /* Function Body */ switch(n__) { case 1: goto L_dscal; } L_dscal: if (*n > 0) { if (*alpha == 0.) { i__1 = (*n - 1) * *incx + 1; i__2 = *incx; for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) { x[ix] = 0.; /* L10: */ } } else if (*alpha == -1.) { i__2 = (*n - 1) * *incx + 1; i__1 = *incx; for (ix = 1; i__1 < 0 ? ix >= i__2 : ix <= i__2; ix += i__1) { x[ix] = -x[ix]; /* L20: */ } } else if (*alpha != 1.) { i__1 = (*n - 1) * *incx + 1; i__2 = *incx; for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) { x[ix] = *alpha * x[ix]; /* L30: */ } } } return 0; /* End of F06EDF. ( DSCAL ) */ } /* f06edf_ */ /* Subroutine */ int f06edf_(n, alpha, x, incx) integer *n; doublereal *alpha, *x; integer *incx; { return f06edf_0_(0, n, alpha, x, incx); } /* Subroutine */ int dscal_(n, alpha, x, incx) integer *n; doublereal *alpha, *x; integer *incx; { return f06edf_0_(1, n, alpha, x, incx); } /* Subroutine */ int f06egf_0_(n__, n, x, incx, y, incy) int n__; integer *n; doublereal *x; integer *incx; doublereal *y; integer *incy; { /* System generated locals */ integer i__1, i__2; /* Local variables */ static doublereal temp; static integer i__, ix, iy; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* .. Entry Points .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* F06EGF performs the operations */ /* temp := x, x := y, y := temp. */ /* Nag Fortran 77 version of the Blas routine DSWAP. */ /* Nag Fortran 77 O( n ) basic linear algebra routine. */ /* -- Written on 26-November-1982. */ /* Sven Hammarling, Nag Central Office. */ /* .. Local Scalars .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --x; --y; /* Function Body */ switch(n__) { case 1: goto L_dswap; } L_dswap: if (*n > 0) { if (*incx == *incy && *incy > 0) { i__1 = (*n - 1) * *incy + 1; i__2 = *incy; for (iy = 1; i__2 < 0 ? iy >= i__1 : iy <= i__1; iy += i__2) { temp = x[iy]; x[iy] = y[iy]; y[iy] = temp; /* L10: */ } } else { if (*incx >= 0) { ix = 1; } else { ix = 1 - (*n - 1) * *incx; } if (*incy > 0) { i__2 = (*n - 1) * *incy + 1; i__1 = *incy; for (iy = 1; i__1 < 0 ? iy >= i__2 : iy <= i__2; iy += i__1) { temp = x[ix]; x[ix] = y[iy]; y[iy] = temp; ix += *incx; /* L20: */ } } else { iy = 1 - (*n - 1) * *incy; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { temp = x[ix]; x[ix] = y[iy]; y[iy] = temp; iy += *incy; ix += *incx; /* L30: */ } } } } return 0; /* End of F06EGF. ( DSWAP ) */ } /* f06egf_ */ /* Subroutine */ int f06egf_(n, x, incx, y, incy) integer *n; doublereal *x; integer *incx; doublereal *y; integer *incy; { return f06egf_0_(0, n, x, incx, y, incy); } /* Subroutine */ int dswap_(n, x, incx, y, incy) integer *n; doublereal *x; integer *incx; doublereal *y; integer *incy; { return f06egf_0_(1, n, x, incx, y, incy); } doublereal f06ejf_0_(n__, n, x, incx) int n__; integer *n; doublereal *x; integer *incx; { /* System generated locals */ doublereal ret_val; /* Local variables */ static doublereal norm; extern doublereal f06bmf_(); extern /* Subroutine */ int f06fjf_(); static doublereal scale, ssq; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* .. Entry Points .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* F06EJF returns the euclidean norm of a vector via the function */ /* name, so that */ /* F06EJF := sqrt( x'*x ) */ /* Nag Fortran 77 version of the Blas routine DNRM2. */ /* Nag Fortran 77 O( n ) basic linear algebra routine. */ /* -- Written on 25-October-1982. */ /* Sven Hammarling, Nag Central Office. */ /* .. Parameters .. */ /* .. Local Scalars .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --x; /* Function Body */ switch(n__) { case 1: goto L_dnrm2; } L_dnrm2: if (*n < 1) { norm = 0.; } else if (*n == 1) { norm = abs(x[1]); } else { scale = 0.; ssq = 1.; f06fjf_(n, &x[1], incx, &scale, &ssq); norm = f06bmf_(&scale, &ssq); } ret_val = norm; return ret_val; /* End of F06EJF. ( DNRM2 ) */ } /* f06ejf_ */ doublereal f06ejf_(n, x, incx) integer *n; doublereal *x; integer *incx; { return f06ejf_0_(0, n, x, incx); } doublereal dnrm2_(n, x, incx) integer *n; doublereal *x; integer *incx; { return f06ejf_0_(1, n, x, incx); } /* Subroutine */ int f06fjf_(n, x, incx, scale, sumsq) integer *n; doublereal *x; integer *incx; doublereal *scale, *sumsq; { /* System generated locals */ integer i__1, i__2; doublereal d__1; /* Local variables */ static doublereal absxi; static integer ix; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* F06FJF returns the values scl and smsq such that */ /* ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq, */ /* where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is assumed */ /* to be at least unity and the value of smsq will then satisfy */ /* 1.0 .le. smsq .le. ( sumsq + n ) . */ /* scale is assumed to be non-negative and scl returns the value */ /* scl = max( scale, abs( x( i ) ) ) . */ /* scale and sumsq must be supplied in SCALE and SUMSQ respectively. */ /* scl and smsq are overwritten on SCALE and SUMSQ respectively. */ /* The routine makes only one pass through the vector X. */ /* Nag Fortran 77 O( n ) basic linear algebra routine. */ /* -- Written on 22-October-1982. */ /* Sven Hammarling, Nag Central Office. */ /* .. Parameters .. */ /* .. Local Scalars .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --x; /* Function Body */ if (*n > 0) { i__1 = (*n - 1) * *incx + 1; i__2 = *incx; for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) { if (x[ix] != 0.) { absxi = (d__1 = x[ix], abs(d__1)); if (*scale < absxi) { /* Computing 2nd power */ d__1 = *scale / absxi; *sumsq = *sumsq * (d__1 * d__1) + 1; *scale = absxi; } else { /* Computing 2nd power */ d__1 = absxi / *scale; *sumsq += d__1 * d__1; } } /* L10: */ } } return 0; /* End of F06FJF. ( SSSQ ) */ } /* f06fjf_ */ /* Subroutine */ int f06frf_(n, alpha, x, incx, tol, zeta) integer *n; doublereal *alpha, *x; integer *incx; doublereal *tol, *zeta; { /* Initialized data */ static logical first = TRUE_; /* System generated locals */ doublereal d__1; /* Builtin functions */ double d_sign(), sqrt(); /* Local variables */ static doublereal beta; extern /* Subroutine */ int f06fjf_(), dscal_(); static doublereal scale; extern doublereal x02ajf_(); static doublereal eps, ssq; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* F06FRF generates details of a generalized Householder reflection such */ /* that */ /* P*( alpha ) = ( beta ), P'*P = I. */ /* ( x ) ( 0 ) */ /* P is given in the form */ /* P = I - ( zeta )*( zeta z' ), */ /* ( z ) */ /* where z is an n element vector and zeta is a scalar that satisfies */ /* 1.0 .le. zeta .le. sqrt( 2.0 ). */ /* zeta is returned in ZETA unless x is such that */ /* max( abs( x( i ) ) ) .le. max( eps*abs( alpha ), tol ) */ /* where eps is the relative machine precision and tol is the user */ /* supplied value TOL, in which case ZETA is returned as 0.0 and P can */ /* be taken to be the unit matrix. */ /* beta is overwritten on alpha and z is overwritten on x. */ /* the routine may be called with n = 0 and advantage is taken of the */ /* case where n = 1. */ /* Nag Fortran 77 O( n ) basic linear algebra routine. */ /* -- Written on 30-August-1984. */ /* Sven Hammarling, Nag Central Office. */ /* This version dated 28-September-1984. */ /* .. Parameters .. */ /* .. Local Scalars .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Save statement .. */ /* .. Data statements .. */ /* Parameter adjustments */ --x; /* Function Body */ /* .. */ /* .. Executable Statements .. */ if (*n < 1) { *zeta = 0.; } else if (*n == 1 && x[1] == 0.) { *zeta = 0.; } else { if (first) { first = FALSE_; eps = x02ajf_(); } /* Treat case where P is a 2 by 2 matrix specially. */ if (*n == 1) { /* Deal with cases where ALPHA = zero and */ /* abs( X( 1 ) ) .le. max( EPS*abs( ALPHA ), TOL ) firs t. */ if (*alpha == 0.) { *zeta = 1.; *alpha = abs(x[1]); x[1] = -d_sign(&c_b15, &x[1]); } else /* if(complicated condition) */ { /* Computing MAX */ d__1 = eps * abs(*alpha); if (abs(x[1]) <= max(d__1,*tol)) { *zeta = 0.; } else { if (abs(*alpha) >= abs(x[1])) { /* Computing 2nd power */ d__1 = x[1] / *alpha; beta = abs(*alpha) * sqrt(d__1 * d__1 + 1); } else { /* Computing 2nd power */ d__1 = *alpha / x[1]; beta = abs(x[1]) * sqrt(d__1 * d__1 + 1); } *zeta = sqrt((abs(*alpha) + beta) / beta); if (*alpha >= 0.) { beta = -beta; } x[1] = -x[1] / (*zeta * beta); *alpha = beta; } } } else { /* Now P is larger than 2 by 2. */ ssq = 1.; scale = 0.; f06fjf_(n, &x[1], incx, &scale, &ssq); /* Treat cases where SCALE = zero, */ /* SCALE .le. max( EPS*abs( ALPHA ), TOL ) and */ /* ALPHA = zero specially. */ /* Note that SCALE = max( abs( X( i ) ) ). */ /* Computing MAX */ d__1 = eps * abs(*alpha); if (scale == 0. || scale <= max(d__1,*tol)) { *zeta = 0.; } else if (*alpha == 0.) { *zeta = 1.; *alpha = scale * sqrt(ssq); d__1 = -1 / *alpha; dscal_(n, &d__1, &x[1], incx); } else { if (scale < abs(*alpha)) { /* Computing 2nd power */ d__1 = scale / *alpha; beta = abs(*alpha) * sqrt(ssq * (d__1 * d__1) + 1); } else { /* Computing 2nd power */ d__1 = *alpha / scale; beta = scale * sqrt(ssq + d__1 * d__1); } *zeta = sqrt((beta + abs(*alpha)) / beta); if (*alpha > 0.) { beta = -beta; } d__1 = -1 / (*zeta * beta); dscal_(n, &d__1, &x[1], incx); *alpha = beta; } } } return 0; /* End of F06FRF. ( SGRFG ) */ } /* f06frf_ */ /* Subroutine */ int f06paf_0_(n__, trans, m, n, alpha, a, lda, x, incx, beta, y, incy, trans_len) int n__; char *trans; integer *m, *n; doublereal *alpha, *a; integer *lda; doublereal *x; integer *incx; doublereal *beta, *y; integer *incy; ftnlen trans_len; { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2; /* Local variables */ static integer info; static doublereal temp; static integer lenx, leny, i__, j; extern /* Subroutine */ int f06aaz_(); static integer ix, iy, jx, jy, kx, ky; /* MARK 13 RE-ISSUE. NAG COPYRIGHT 1988. */ /* .. Entry Points .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* DGEMV performs one of the matrix-vector operations */ /* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, */ /* where alpha and beta are scalars, x and y are vectors and A is an */ /* m by n matrix. */ /* Parameters */ /* ========== */ /* TRANS - CHARACTER*1. */ /* On entry, TRANS specifies the operation to be performed as */ /* follows: */ /* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */ /* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */ /* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. */ /* Unchanged on exit. */ /* M - INTEGER. */ /* On entry, M specifies the number of rows of the matrix A. */ /* M must be at least zero. */ /* Unchanged on exit. */ /* N - INTEGER. */ /* On entry, N specifies the number of columns of the matrix A. */ /* N must be at least zero. */ /* Unchanged on exit. */ /* ALPHA - DOUBLE PRECISION. */ /* On entry, ALPHA specifies the scalar alpha. */ /* Unchanged on exit. */ /* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */ /* Before entry, the leading m by n part of the array A must */ /* contain the matrix of coefficients. */ /* Unchanged on exit. */ /* LDA - INTEGER. */ /* On entry, LDA specifies the first dimension of A as declared */ /* in the calling (sub) program. LDA must be at least */ /* max( 1, m ). */ /* Unchanged on exit. */ /* X - DOUBLE PRECISION array of DIMENSION at least */ /* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */ /* and at least */ /* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */ /* Before entry, the incremented array X must contain the */ /* vector x. */ /* Unchanged on exit. */ /* INCX - INTEGER. */ /* On entry, INCX specifies the increment for the elements of */ /* X. INCX must not be zero. */ /* Unchanged on exit. */ /* BETA - DOUBLE PRECISION. */ /* On entry, BETA specifies the scalar beta. When BETA is */ /* supplied as zero then Y need not be set on input. */ /* Unchanged on exit. */ /* Y - DOUBLE PRECISION array of DIMENSION at least */ /* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */ /* and at least */ /* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */ /* Before entry with BETA non-zero, the incremented array Y */ /* must contain the vector y. On exit, Y is overwritten by the */ /* updated vector y. */ /* INCY - INTEGER. */ /* On entry, INCY specifies the increment for the elements of */ /* Y. INCY must not be zero. */ /* Unchanged on exit. */ /* Level 2 Blas routine. */ /* -- Written on 22-October-1986. */ /* Jack Dongarra, Argonne National Lab. */ /* Jeremy Du Croz, Nag Central Office. */ /* Sven Hammarling, Nag Central Office. */ /* Richard Hanson, Sandia National Labs. */ /* .. Parameters .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Test the input parameters. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; --x; --y; /* Function Body */ switch(n__) { case 1: goto L_dgemv; } L_dgemv: info = 0; if (! (*(unsigned char *)trans == 'N' || *(unsigned char *)trans == 'n') && ! (*(unsigned char *)trans == 'T' || *(unsigned char *)trans == 't') && ! (*(unsigned char *)trans == 'C' || *(unsigned char *) trans == 'c')) { info = 1; } else if (*m < 0) { info = 2; } else if (*n < 0) { info = 3; } else if (*lda < max(1,*m)) { info = 6; } else if (*incx == 0) { info = 8; } else if (*incy == 0) { info = 11; } if (info != 0) { f06aaz_("F06PAF/DGEMV ", &info, 13L); return 0; } /* Quick return if possible. */ if (*m == 0 || *n == 0 || *alpha == 0. && *beta == 1.) { return 0; } /* Set LENX and LENY, the lengths of the vectors x and y, and set */ /* up the start points in X and Y. */ if (*(unsigned char *)trans == 'N' || *(unsigned char *)trans == 'n') { lenx = *n; leny = *m; } else { lenx = *m; leny = *n; } if (*incx > 0) { kx = 1; } else { kx = 1 - (lenx - 1) * *incx; } if (*incy > 0) { ky = 1; } else { ky = 1 - (leny - 1) * *incy; } /* Start the operations. In this version the elements of A are */ /* accessed sequentially with one pass through A. */ /* First form y := beta*y. */ if (*beta != 1.) { if (*incy == 1) { if (*beta == 0.) { i__1 = leny; for (i__ = 1; i__ <= i__1; ++i__) { y[i__] = 0.; /* L10: */ } } else { i__1 = leny; for (i__ = 1; i__ <= i__1; ++i__) { y[i__] = *beta * y[i__]; /* L20: */ } } } else { iy = ky; if (*beta == 0.) { i__1 = leny; for (i__ = 1; i__ <= i__1; ++i__) { y[iy] = 0.; iy += *incy; /* L30: */ } } else { i__1 = leny; for (i__ = 1; i__ <= i__1; ++i__) { y[iy] = *beta * y[iy]; iy += *incy; /* L40: */ } } } } if (*alpha == 0.) { return 0; } if (*(unsigned char *)trans == 'N' || *(unsigned char *)trans == 'n') { /* Form y := alpha*A*x + y. */ jx = kx; if (*incy == 1) { i__1 = *n; for (j = 1; j <= i__1; ++j) { if (x[jx] != 0.) { temp = *alpha * x[jx]; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { y[i__] += temp * a[i__ + j * a_dim1]; /* L50: */ } } jx += *incx; /* L60: */ } } else { i__1 = *n; for (j = 1; j <= i__1; ++j) { if (x[jx] != 0.) { temp = *alpha * x[jx]; iy = ky; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { y[iy] += temp * a[i__ + j * a_dim1]; iy += *incy; /* L70: */ } } jx += *incx; /* L80: */ } } } else { /* Form y := alpha*A'*x + y. */ jy = ky; if (*incx == 1) { i__1 = *n; for (j = 1; j <= i__1; ++j) { temp = 0.; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { temp += a[i__ + j * a_dim1] * x[i__]; /* L90: */ } y[jy] += *alpha * temp; jy += *incy; /* L100: */ } } else { i__1 = *n; for (j = 1; j <= i__1; ++j) { temp = 0.; ix = kx; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { temp += a[i__ + j * a_dim1] * x[ix]; ix += *incx; /* L110: */ } y[jy] += *alpha * temp; jy += *incy; /* L120: */ } } } return 0; /* End of F06PAF (DGEMV ). */ } /* f06paf_ */ /* Subroutine */ int f06paf_(trans, m, n, alpha, a, lda, x, incx, beta, y, incy, trans_len) char *trans; integer *m, *n; doublereal *alpha, *a; integer *lda; doublereal *x; integer *incx; doublereal *beta, *y; integer *incy; ftnlen trans_len; { return f06paf_0_(0, trans, m, n, alpha, a, lda, x, incx, beta, y, incy, trans_len); } /* Subroutine */ int dgemv_(trans, m, n, alpha, a, lda, x, incx, beta, y, incy, trans_len) char *trans; integer *m, *n; doublereal *alpha, *a; integer *lda; doublereal *x; integer *incx; doublereal *beta, *y; integer *incy; ftnlen trans_len; { return f06paf_0_(1, trans, m, n, alpha, a, lda, x, incx, beta, y, incy, trans_len); } /* Subroutine */ int f06pjf_0_(n__, uplo, trans, diag, n, a, lda, x, incx, uplo_len, trans_len, diag_len) int n__; char *uplo, *trans, *diag; integer *n; doublereal *a; integer *lda; doublereal *x; integer *incx; ftnlen uplo_len; ftnlen trans_len; ftnlen diag_len; { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2; /* Local variables */ static integer info; static doublereal temp; static integer i__, j; extern /* Subroutine */ int f06aaz_(); static integer ix, jx, kx; static logical nounit; /* MARK 13 RE-ISSUE. NAG COPYRIGHT 1988. */ /* .. Entry Points .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* DTRSV solves one of the systems of equations */ /* A*x = b, or A'*x = b, */ /* where b and x are n element vectors and A is an n by n unit, or */ /* non-unit, upper or lower triangular matrix. */ /* No test for singularity or near-singularity is included in this */ /* routine. Such tests must be performed before calling this routine. */ /* Parameters */ /* ========== */ /* UPLO - CHARACTER*1. */ /* On entry, UPLO specifies whether the matrix is an upper or */ /* lower triangular matrix as follows: */ /* UPLO = 'U' or 'u' A is an upper triangular matrix. */ /* UPLO = 'L' or 'l' A is a lower triangular matrix. */ /* Unchanged on exit. */ /* TRANS - CHARACTER*1. */ /* On entry, TRANS specifies the equations to be solved as */ /* follows: */ /* TRANS = 'N' or 'n' A*x = b. */ /* TRANS = 'T' or 't' A'*x = b. */ /* TRANS = 'C' or 'c' A'*x = b. */ /* Unchanged on exit. */ /* DIAG - CHARACTER*1. */ /* On entry, DIAG specifies whether or not A is unit */ /* triangular as follows: */ /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */ /* DIAG = 'N' or 'n' A is not assumed to be unit */ /* triangular. */ /* Unchanged on exit. */ /* N - INTEGER. */ /* On entry, N specifies the order of the matrix A. */ /* N must be at least zero. */ /* Unchanged on exit. */ /* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */ /* Before entry with UPLO = 'U' or 'u', the leading n by n */ /* upper triangular part of the array A must contain the upper */ /* triangular matrix and the strictly lower triangular part of */ /* A is not referenced. */ /* Before entry with UPLO = 'L' or 'l', the leading n by n */ /* lower triangular part of the array A must contain the lower */ /* triangular matrix and the strictly upper triangular part of */ /* A is not referenced. */ /* Note that when DIAG = 'U' or 'u', the diagonal elements of */ /* A are not referenced either, but are assumed to be unity. */ /* Unchanged on exit. */ /* LDA - INTEGER. */ /* On entry, LDA specifies the first dimension of A as declared */ /* in the calling (sub) program. LDA must be at least */ /* max( 1, n ). */ /* Unchanged on exit. */ /* X - DOUBLE PRECISION array of dimension at least */ /* ( 1 + ( n - 1 )*abs( INCX ) ). */ /* Before entry, the incremented array X must contain the n */ /* element right-hand side vector b. On exit, X is overwritten */ /* with the solution vector x. */ /* INCX - INTEGER. */ /* On entry, INCX specifies the increment for the elements of */ /* X. INCX must not be zero. */ /* Unchanged on exit. */ /* Level 2 Blas routine. */ /* -- Written on 22-October-1986. */ /* Jack Dongarra, Argonne National Lab. */ /* Jeremy Du Croz, Nag Central Office. */ /* Sven Hammarling, Nag Central Office. */ /* Richard Hanson, Sandia National Labs. */ /* .. Parameters .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Test the input parameters. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; --x; /* Function Body */ switch(n__) { case 1: goto L_dtrsv; } L_dtrsv: info = 0; if (! (*(unsigned char *)uplo == 'U' || *(unsigned char *)uplo == 'u') && ! (*(unsigned char *)uplo == 'L' || *(unsigned char *)uplo == 'l') ) { info = 1; } else if (! (*(unsigned char *)trans == 'N' || *(unsigned char *)trans == 'n') && ! (*(unsigned char *)trans == 'T' || *(unsigned char *) trans == 't') && ! (*(unsigned char *)trans == 'C' || *(unsigned char *)trans == 'c')) { info = 2; } else if (! (*(unsigned char *)diag == 'U' || *(unsigned char *)diag == 'u') && ! (*(unsigned char *)diag == 'N' || *(unsigned char *) diag == 'n')) { info = 3; } else if (*n < 0) { info = 4; } else if (*lda < max(1,*n)) { info = 6; } else if (*incx == 0) { info = 8; } if (info != 0) { f06aaz_("F06PJF/DTRSV ", &info, 13L); return 0; } /* Quick return if possible. */ if (*n == 0) { return 0; } nounit = *(unsigned char *)diag == 'N' || *(unsigned char *)diag == 'n'; /* Set up the start point in X if the increment is not unity. This */ /* will be ( N - 1 )*INCX too small for descending loops. */ if (*incx <= 0) { kx = 1 - (*n - 1) * *incx; } else if (*incx != 1) { kx = 1; } /* Start the operations. In this version the elements of A are */ /* accessed sequentially with one pass through A. */ if (*(unsigned char *)trans == 'N' || *(unsigned char *)trans == 'n') { /* Form x := inv( A )*x. */ if (*(unsigned char *)uplo == 'U' || *(unsigned char *)uplo == 'u') { if (*incx == 1) { for (j = *n; j >= 1; --j) { if (x[j] != 0.) { if (nounit) { x[j] /= a[j + j * a_dim1]; } temp = x[j]; for (i__ = j - 1; i__ >= 1; --i__) { x[i__] -= temp * a[i__ + j * a_dim1]; /* L10: */ } } /* L20: */ } } else { jx = kx + (*n - 1) * *incx; for (j = *n; j >= 1; --j) { if (x[jx] != 0.) { if (nounit) { x[jx] /= a[j + j * a_dim1]; } temp = x[jx]; ix = jx; for (i__ = j - 1; i__ >= 1; --i__) { ix -= *incx; x[ix] -= temp * a[i__ + j * a_dim1]; /* L30: */ } } jx -= *incx; /* L40: */ } } } else { if (*incx == 1) { i__1 = *n; for (j = 1; j <= i__1; ++j) { if (x[j] != 0.) { if (nounit) { x[j] /= a[j + j * a_dim1]; } temp = x[j]; i__2 = *n; for (i__ = j + 1; i__ <= i__2; ++i__) { x[i__] -= temp * a[i__ + j * a_dim1]; /* L50: */ } } /* L60: */ } } else { jx = kx; i__1 = *n; for (j = 1; j <= i__1; ++j) { if (x[jx] != 0.) { if (nounit) { x[jx] /= a[j + j * a_dim1]; } temp = x[jx]; ix = jx; i__2 = *n; for (i__ = j + 1; i__ <= i__2; ++i__) { ix += *incx; x[ix] -= temp * a[i__ + j * a_dim1]; /* L70: */ } } jx += *incx; /* L80: */ } } } } else { /* Form x := inv( A' )*x. */ if (*(unsigned char *)uplo == 'U' || *(unsigned char *)uplo == 'u') { if (*incx == 1) { i__1 = *n; for (j = 1; j <= i__1; ++j) { temp = x[j]; i__2 = j - 1; for (i__ = 1; i__ <= i__2; ++i__) { temp -= a[i__ + j * a_dim1] * x[i__]; /* L90: */ } if (nounit) { temp /= a[j + j * a_dim1]; } x[j] = temp; /* L100: */ } } else { jx = kx; i__1 = *n; for (j = 1; j <= i__1; ++j) { temp = x[jx]; ix = kx; i__2 = j - 1; for (i__ = 1; i__ <= i__2; ++i__) { temp -= a[i__ + j * a_dim1] * x[ix]; ix += *incx; /* L110: */ } if (nounit) { temp /= a[j + j * a_dim1]; } x[jx] = temp; jx += *incx; /* L120: */ } } } else { if (*incx == 1) { for (j = *n; j >= 1; --j) { temp = x[j]; i__1 = j + 1; for (i__ = *n; i__ >= i__1; --i__) { temp -= a[i__ + j * a_dim1] * x[i__]; /* L130: */ } if (nounit) { temp /= a[j + j * a_dim1]; } x[j] = temp; /* L140: */ } } else { kx += (*n - 1) * *incx; jx = kx; for (j = *n; j >= 1; --j) { temp = x[jx]; ix = kx; i__1 = j + 1; for (i__ = *n; i__ >= i__1; --i__) { temp -= a[i__ + j * a_dim1] * x[ix]; ix -= *incx; /* L150: */ } if (nounit) { temp /= a[j + j * a_dim1]; } x[jx] = temp; jx -= *incx; /* L160: */ } } } } return 0; /* End of F06PJF (DTRSV ). */ } /* f06pjf_ */ /* Subroutine */ int f06pjf_(uplo, trans, diag, n, a, lda, x, incx, uplo_len, trans_len, diag_len) char *uplo, *trans, *diag; integer *n; doublereal *a; integer *lda; doublereal *x; integer *incx; ftnlen uplo_len; ftnlen trans_len; ftnlen diag_len; { return f06pjf_0_(0, uplo, trans, diag, n, a, lda, x, incx, uplo_len, trans_len, diag_len); } /* Subroutine */ int dtrsv_(uplo, trans, diag, n, a, lda, x, incx, uplo_len, trans_len, diag_len) char *uplo, *trans, *diag; integer *n; doublereal *a; integer *lda; doublereal *x; integer *incx; ftnlen uplo_len; ftnlen trans_len; ftnlen diag_len; { return f06pjf_0_(1, uplo, trans, diag, n, a, lda, x, incx, uplo_len, trans_len, diag_len); } /* Subroutine */ int f06pmf_0_(n__, m, n, alpha, x, incx, y, incy, a, lda) int n__; integer *m, *n; doublereal *alpha, *x; integer *incx; doublereal *y; integer *incy; doublereal *a; integer *lda; { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2; /* Local variables */ static integer info; static doublereal temp; static integer i__, j; extern /* Subroutine */ int f06aaz_(); static integer ix, jy, kx; /* MARK 13 RE-ISSUE. NAG COPYRIGHT 1988. */ /* .. Entry Points .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* DGER performs the rank 1 operation */ /* A := alpha*x*y' + A, */ /* where alpha is a scalar, x is an m element vector, y is an n element */ /* vector and A is an m by n matrix. */ /* Parameters */ /* ========== */ /* M - INTEGER. */ /* On entry, M specifies the number of rows of the matrix A. */ /* M must be at least zero. */ /* Unchanged on exit. */ /* N - INTEGER. */ /* On entry, N specifies the number of columns of the matrix A. */ /* N must be at least zero. */ /* Unchanged on exit. */ /* ALPHA - DOUBLE PRECISION. */ /* On entry, ALPHA specifies the scalar alpha. */ /* Unchanged on exit. */ /* X - DOUBLE PRECISION array of dimension at least */ /* ( 1 + ( m - 1 )*abs( INCX ) ). */ /* Before entry, the incremented array X must contain the m */ /* element vector x. */ /* Unchanged on exit. */ /* INCX - INTEGER. */ /* On entry, INCX specifies the increment for the elements of */ /* X. INCX must not be zero. */ /* Unchanged on exit. */ /* Y - DOUBLE PRECISION array of dimension at least */ /* ( 1 + ( n - 1 )*abs( INCY ) ). */ /* Before entry, the incremented array Y must contain the n */ /* element vector y. */ /* Unchanged on exit. */ /* INCY - INTEGER. */ /* On entry, INCY specifies the increment for the elements of */ /* Y. INCY must not be zero. */ /* Unchanged on exit. */ /* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */ /* Before entry, the leading m by n part of the array A must */ /* contain the matrix of coefficients. On exit, A is */ /* overwritten by the updated matrix. */ /* LDA - INTEGER. */ /* On entry, LDA specifies the first dimension of A as declared */ /* in the calling (sub) program. LDA must be at least */ /* max( 1, m ). */ /* Unchanged on exit. */ /* Level 2 Blas routine. */ /* -- Written on 22-October-1986. */ /* Jack Dongarra, Argonne National Lab. */ /* Jeremy Du Croz, Nag Central Office. */ /* Sven Hammarling, Nag Central Office. */ /* Richard Hanson, Sandia National Labs. */ /* .. Parameters .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Executable Statements .. */ /* Test the input parameters. */ /* Parameter adjustments */ --x; --y; a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ switch(n__) { case 1: goto L_dger; } L_dger: info = 0; if (*m < 0) { info = 1; } else if (*n < 0) { info = 2; } else if (*incx == 0) { info = 5; } else if (*incy == 0) { info = 7; } else if (*lda < max(1,*m)) { info = 9; } if (info != 0) { f06aaz_("F06PMF/DGER ", &info, 13L); return 0; } /* Quick return if possible. */ if (*m == 0 || *n == 0 || *alpha == 0.) { return 0; } /* Start the operations. In this version the elements of A are */ /* accessed sequentially with one pass through A. */ if (*incy > 0) { jy = 1; } else { jy = 1 - (*n - 1) * *incy; } if (*incx == 1) { i__1 = *n; for (j = 1; j <= i__1; ++j) { if (y[jy] != 0.) { temp = *alpha * y[jy]; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { a[i__ + j * a_dim1] += x[i__] * temp; /* L10: */ } } jy += *incy; /* L20: */ } } else { if (*incx > 0) { kx = 1; } else { kx = 1 - (*m - 1) * *incx; } i__1 = *n; for (j = 1; j <= i__1; ++j) { if (y[jy] != 0.) { temp = *alpha * y[jy]; ix = kx; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { a[i__ + j * a_dim1] += x[ix] * temp; ix += *incx; /* L30: */ } } jy += *incy; /* L40: */ } } return 0; /* End of F06PMF (DGER ). */ } /* f06pmf_ */ /* Subroutine */ int f06pmf_(m, n, alpha, x, incx, y, incy, a, lda) integer *m, *n; doublereal *alpha, *x; integer *incx; doublereal *y; integer *incy; doublereal *a; integer *lda; { return f06pmf_0_(0, m, n, alpha, x, incx, y, incy, a, lda); } /* Subroutine */ int dger_(m, n, alpha, x, incx, y, incy, a, lda) integer *m, *n; doublereal *alpha, *x; integer *incx; doublereal *y; integer *incy; doublereal *a; integer *lda; { return f06pmf_0_(1, m, n, alpha, x, incx, y, incy, a, lda); } /* Subroutine */ int f06yaf_0_(n__, transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c__, ldc, transa_len, transb_len) int n__; char *transa, *transb; integer *m, *n, *k; doublereal *alpha, *a; integer *lda; doublereal *b; integer *ldb; doublereal *beta, *c__; integer *ldc; ftnlen transa_len; ftnlen transb_len; { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, i__3; /* Local variables */ static integer info; static logical nota, notb; static doublereal temp; static integer i__, j, l; extern /* Subroutine */ int f06aaz_(); static integer ncola, nrowa, nrowb; /* MARK 14 RELEASE. NAG COPYRIGHT 1989. */ /* Jack Dongarra, Argonne National Laboratory. */ /* Iain Duff, AERE Harwell. */ /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */ /* Sven Hammarling, Numerical Algorithms Group Ltd. */ /* .. Entry Points .. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Set NOTA and NOTB as true if A and B respectively are not */ /* transposed and set NROWA, NCOLA and NROWB as the number of rows */ /* and columns of A and the number of rows of B respectively. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; b_dim1 = *ldb; b_offset = b_dim1 + 1; b -= b_offset; c_dim1 = *ldc; c_offset = c_dim1 + 1; c__ -= c_offset; /* Function Body */ switch(n__) { case 1: goto L_dgemm; } L_dgemm: nota = *(unsigned char *)transa == 'N' || *(unsigned char *)transa == 'n'; notb = *(unsigned char *)transb == 'N' || *(unsigned char *)transb == 'n'; if (nota) { nrowa = *m; ncola = *k; } else { nrowa = *k; ncola = *m; } if (notb) { nrowb = *k; } else { nrowb = *n; } /* Test the input parameters. */ info = 0; if (! nota && ! (*(unsigned char *)transa == 'C' || *(unsigned char *) transa == 'c') && ! (*(unsigned char *)transa == 'T' || *( unsigned char *)transa == 't')) { info = 1; } else if (! notb && ! (*(unsigned char *)transb == 'C' || *(unsigned char *)transb == 'c') && ! (*(unsigned char *)transb == 'T' || *( unsigned char *)transb == 't')) { info = 2; } else if (*m < 0) { info = 3; } else if (*n < 0) { info = 4; } else if (*k < 0) { info = 5; } else if (*lda < max(1,nrowa)) { info = 8; } else if (*ldb < max(1,nrowb)) { info = 10; } else if (*ldc < max(1,*m)) { info = 13; } if (info != 0) { f06aaz_("F06YAF/DGEMM ", &info, 13L); return 0; } /* Quick return if possible. */ if (*m == 0 || *n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.) { return 0; } /* And if alpha.eq.zero. */ if (*alpha == 0.) { if (*beta == 0.) { i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] = 0.; /* L20: */ } /* L40: */ } } else { i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]; /* L60: */ } /* L80: */ } } return 0; } /* Start the operations. */ if (notb) { if (nota) { /* Form C := alpha*A*B + beta*C. */ i__1 = *n; for (j = 1; j <= i__1; ++j) { if (*beta == 0.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] = 0.; /* L100: */ } } else if (*beta != 1.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]; /* L120: */ } } i__2 = *k; for (l = 1; l <= i__2; ++l) { if (b[l + j * b_dim1] != 0.) { temp = *alpha * b[l + j * b_dim1]; i__3 = *m; for (i__ = 1; i__ <= i__3; ++i__) { c__[i__ + j * c_dim1] += temp * a[i__ + l * a_dim1]; /* L140: */ } } /* L160: */ } /* L180: */ } } else { /* Form C := alpha*A'*B + beta*C */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { temp = 0.; i__3 = *k; for (l = 1; l <= i__3; ++l) { temp += a[l + i__ * a_dim1] * b[l + j * b_dim1]; /* L200: */ } if (*beta == 0.) { c__[i__ + j * c_dim1] = *alpha * temp; } else { c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[ i__ + j * c_dim1]; } /* L220: */ } /* L240: */ } } } else { if (nota) { /* Form C := alpha*A*B' + beta*C */ i__1 = *n; for (j = 1; j <= i__1; ++j) { if (*beta == 0.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] = 0.; /* L260: */ } } else if (*beta != 1.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1]; /* L280: */ } } i__2 = *k; for (l = 1; l <= i__2; ++l) { if (b[j + l * b_dim1] != 0.) { temp = *alpha * b[j + l * b_dim1]; i__3 = *m; for (i__ = 1; i__ <= i__3; ++i__) { c__[i__ + j * c_dim1] += temp * a[i__ + l * a_dim1]; /* L300: */ } } /* L320: */ } /* L340: */ } } else { /* Form C := alpha*A'*B' + beta*C */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { temp = 0.; i__3 = *k; for (l = 1; l <= i__3; ++l) { temp += a[l + i__ * a_dim1] * b[j + l * b_dim1]; /* L360: */ } if (*beta == 0.) { c__[i__ + j * c_dim1] = *alpha * temp; } else { c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[ i__ + j * c_dim1]; } /* L380: */ } /* L400: */ } } } return 0; /* End of F06YAF (DGEMM ). */ } /* f06yaf_ */ /* Subroutine */ int f06yaf_(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c__, ldc, transa_len, transb_len) char *transa, *transb; integer *m, *n, *k; doublereal *alpha, *a; integer *lda; doublereal *b; integer *ldb; doublereal *beta, *c__; integer *ldc; ftnlen transa_len; ftnlen transb_len; { return f06yaf_0_(0, transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c__, ldc, transa_len, transb_len); } /* Subroutine */ int dgemm_(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c__, ldc, transa_len, transb_len) char *transa, *transb; integer *m, *n, *k; doublereal *alpha, *a; integer *lda; doublereal *b; integer *ldb; doublereal *beta, *c__; integer *ldc; ftnlen transa_len; ftnlen transb_len; { return f06yaf_0_(1, transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c__, ldc, transa_len, transb_len); } /* Subroutine */ int f06yjf_0_(n__, side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb, side_len, uplo_len, transa_len, diag_len) int n__; char *side, *uplo, *transa, *diag; integer *m, *n; doublereal *alpha, *a; integer *lda; doublereal *b; integer *ldb; ftnlen side_len; ftnlen uplo_len; ftnlen transa_len; ftnlen diag_len; { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3; /* Local variables */ static integer info; static doublereal temp; static integer i__, j, k; extern /* Subroutine */ int f06aaz_(); static logical lside; static integer nrowa; static logical upper, nounit; /* MARK 14 RELEASE. NAG COPYRIGHT 1989. */ /* Purpose */ /* ======= */ /* DTRSM solves one of the matrix equations */ /* op( A )*X = alpha*B, or X*op( A ) = alpha*B, */ /* where alpha is a scalar, X and B are m by n matrices, A is a unit, or */ /* non-unit, upper or lower triangular matrix and op( A ) is one of */ /* op( A ) = A or op( A ) = A'. */ /* The matrix X is overwritten on B. */ /* Parameters */ /* ========== */ /* SIDE - CHARACTER*1. */ /* On entry, SIDE specifies whether op( A ) appears on the left */ /* or right of X as follows: */ /* SIDE = 'L' or 'l' op( A )*X = alpha*B. */ /* SIDE = 'R' or 'r' X*op( A ) = alpha*B. */ /* Unchanged on exit. */ /* UPLO - CHARACTER*1. */ /* On entry, UPLO specifies whether the matrix A is an upper or */ /* lower triangular matrix as follows: */ /* UPLO = 'U' or 'u' A is an upper triangular matrix. */ /* UPLO = 'L' or 'l' A is a lower triangular matrix. */ /* Unchanged on exit. */ /* TRANSA - CHARACTER*1. */ /* On entry, TRANSA specifies the form of op( A ) to be used in */ /* the matrix multiplication as follows: */ /* TRANSA = 'N' or 'n' op( A ) = A. */ /* TRANSA = 'T' or 't' op( A ) = A'. */ /* TRANSA = 'C' or 'c' op( A ) = A'. */ /* Unchanged on exit. */ /* DIAG - CHARACTER*1. */ /* On entry, DIAG specifies whether or not A is unit triangular */ /* as follows: */ /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */ /* DIAG = 'N' or 'n' A is not assumed to be unit */ /* triangular. */ /* Unchanged on exit. */ /* M - INTEGER. */ /* On entry, M specifies the number of rows of B. M must be at */ /* least zero. */ /* Unchanged on exit. */ /* N - INTEGER. */ /* On entry, N specifies the number of columns of B. N must be */ /* at least zero. */ /* Unchanged on exit. */ /* ALPHA - REAL . */ /* On entry, ALPHA specifies the scalar alpha. When alpha is */ /* zero then A is not referenced and B need not be set before */ /* entry. */ /* Unchanged on exit. */ /* A - REAL array of DIMENSION ( LDA, k ), where k is m */ /* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. */ /* Before entry with UPLO = 'U' or 'u', the leading k by k */ /* upper triangular part of the array A must contain the upper */ /* triangular matrix and the strictly lower triangular part of */ /* A is not referenced. */ /* Before entry with UPLO = 'L' or 'l', the leading k by k */ /* lower triangular part of the array A must contain the lower */ /* triangular matrix and the strictly upper triangular part of */ /* A is not referenced. */ /* Note that when DIAG = 'U' or 'u', the diagonal elements of */ /* A are not referenced either, but are assumed to be unity. */ /* Unchanged on exit. */ /* LDA - INTEGER. */ /* On entry, LDA specifies the first dimension of A as declared */ /* in the calling (sub) program. When SIDE = 'L' or 'l' then */ /* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' */ /* then LDA must be at least max( 1, n ). */ /* Unchanged on exit. */ /* B - REAL array of DIMENSION ( LDB, n ). */ /* Before entry, the leading m by n part of the array B must */ /* contain the right-hand side matrix B, and on exit is */ /* overwritten by the solution matrix X. */ /* LDB - INTEGER. */ /* On entry, LDB specifies the first dimension of B as declared */ /* in the calling (sub) program. LDB must be at least */ /* max( 1, m ). */ /* Unchanged on exit. */ /* Level 3 Blas routine. */ /* -- Written on 8-February-1989. */ /* Jack Dongarra, Argonne National Laboratory. */ /* Iain Duff, AERE Harwell. */ /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */ /* Sven Hammarling, Numerical Algorithms Group Ltd. */ /* .. Entry Points .. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Test the input parameters. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; b_dim1 = *ldb; b_offset = b_dim1 + 1; b -= b_offset; /* Function Body */ switch(n__) { case 1: goto L_dtrsm; } L_dtrsm: lside = *(unsigned char *)side == 'L' || *(unsigned char *)side == 'l'; if (lside) { nrowa = *m; } else { nrowa = *n; } nounit = *(unsigned char *)diag == 'N' || *(unsigned char *)diag == 'n'; upper = *(unsigned char *)uplo == 'U' || *(unsigned char *)uplo == 'u'; info = 0; if (! lside && ! (*(unsigned char *)side == 'R' || *(unsigned char *)side == 'r')) { info = 1; } else if (! upper && ! (*(unsigned char *)uplo == 'L' || *(unsigned char *)uplo == 'l')) { info = 2; } else if (! (*(unsigned char *)transa == 'N' || *(unsigned char *)transa == 'n') && ! (*(unsigned char *)transa == 'T' || *(unsigned char * )transa == 't') && ! (*(unsigned char *)transa == 'C' || *( unsigned char *)transa == 'c')) { info = 3; } else if (! (*(unsigned char *)diag == 'U' || *(unsigned char *)diag == 'u') && ! (*(unsigned char *)diag == 'N' || *(unsigned char *) diag == 'n')) { info = 4; } else if (*m < 0) { info = 5; } else if (*n < 0) { info = 6; } else if (*lda < max(1,nrowa)) { info = 9; } else if (*ldb < max(1,*m)) { info = 11; } if (info != 0) { f06aaz_("F06YJF/DTRSM ", &info, 13L); return 0; } /* Quick return if possible. */ if (*n == 0) { return 0; } /* And when alpha.eq.zero. */ if (*alpha == 0.) { i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + j * b_dim1] = 0.; /* L20: */ } /* L40: */ } return 0; } /* Start the operations. */ if (lside) { if (*(unsigned char *)transa == 'N' || *(unsigned char *)transa == 'n') { /* Form B := alpha*inv( A )*B. */ if (upper) { i__1 = *n; for (j = 1; j <= i__1; ++j) { if (*alpha != 1.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1] ; /* L60: */ } } for (k = *m; k >= 1; --k) { if (b[k + j * b_dim1] != 0.) { if (nounit) { b[k + j * b_dim1] /= a[k + k * a_dim1]; } i__2 = k - 1; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[ i__ + k * a_dim1]; /* L80: */ } } /* L100: */ } /* L120: */ } } else { i__1 = *n; for (j = 1; j <= i__1; ++j) { if (*alpha != 1.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1] ; /* L140: */ } } i__2 = *m; for (k = 1; k <= i__2; ++k) { if (b[k + j * b_dim1] != 0.) { if (nounit) { b[k + j * b_dim1] /= a[k + k * a_dim1]; } i__3 = *m; for (i__ = k + 1; i__ <= i__3; ++i__) { b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[ i__ + k * a_dim1]; /* L160: */ } } /* L180: */ } /* L200: */ } } } else { /* Form B := alpha*inv( A' )*B. */ if (upper) { i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { temp = *alpha * b[i__ + j * b_dim1]; i__3 = i__ - 1; for (k = 1; k <= i__3; ++k) { temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1]; /* L220: */ } if (nounit) { temp /= a[i__ + i__ * a_dim1]; } b[i__ + j * b_dim1] = temp; /* L240: */ } /* L260: */ } } else { i__1 = *n; for (j = 1; j <= i__1; ++j) { for (i__ = *m; i__ >= 1; --i__) { temp = *alpha * b[i__ + j * b_dim1]; i__2 = *m; for (k = i__ + 1; k <= i__2; ++k) { temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1]; /* L280: */ } if (nounit) { temp /= a[i__ + i__ * a_dim1]; } b[i__ + j * b_dim1] = temp; /* L300: */ } /* L320: */ } } } } else { if (*(unsigned char *)transa == 'N' || *(unsigned char *)transa == 'n') { /* Form B := alpha*B*inv( A ). */ if (upper) { i__1 = *n; for (j = 1; j <= i__1; ++j) { if (*alpha != 1.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1] ; /* L340: */ } } i__2 = j - 1; for (k = 1; k <= i__2; ++k) { if (a[k + j * a_dim1] != 0.) { i__3 = *m; for (i__ = 1; i__ <= i__3; ++i__) { b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[ i__ + k * b_dim1]; /* L360: */ } } /* L380: */ } if (nounit) { temp = 1. / a[j + j * a_dim1]; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1]; /* L400: */ } } /* L420: */ } } else { for (j = *n; j >= 1; --j) { if (*alpha != 1.) { i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1] ; /* L440: */ } } i__1 = *n; for (k = j + 1; k <= i__1; ++k) { if (a[k + j * a_dim1] != 0.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[ i__ + k * b_dim1]; /* L460: */ } } /* L480: */ } if (nounit) { temp = 1. / a[j + j * a_dim1]; i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1]; /* L500: */ } } /* L520: */ } } } else { /* Form B := alpha*B*inv( A' ). */ if (upper) { for (k = *n; k >= 1; --k) { if (nounit) { temp = 1. / a[k + k * a_dim1]; i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1]; /* L540: */ } } i__1 = k - 1; for (j = 1; j <= i__1; ++j) { if (a[j + k * a_dim1] != 0.) { temp = a[j + k * a_dim1]; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + j * b_dim1] -= temp * b[i__ + k * b_dim1]; /* L560: */ } } /* L580: */ } if (*alpha != 1.) { i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1] ; /* L600: */ } } /* L620: */ } } else { i__1 = *n; for (k = 1; k <= i__1; ++k) { if (nounit) { temp = 1. / a[k + k * a_dim1]; i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1]; /* L640: */ } } i__2 = *n; for (j = k + 1; j <= i__2; ++j) { if (a[j + k * a_dim1] != 0.) { temp = a[j + k * a_dim1]; i__3 = *m; for (i__ = 1; i__ <= i__3; ++i__) { b[i__ + j * b_dim1] -= temp * b[i__ + k * b_dim1]; /* L660: */ } } /* L680: */ } if (*alpha != 1.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1] ; /* L700: */ } } /* L720: */ } } } } return 0; /* End of F06YJF (DTRSM ). */ } /* f06yjf_ */ /* Subroutine */ int f06yjf_(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb, side_len, uplo_len, transa_len, diag_len) char *side, *uplo, *transa, *diag; integer *m, *n; doublereal *alpha, *a; integer *lda; doublereal *b; integer *ldb; ftnlen side_len; ftnlen uplo_len; ftnlen transa_len; ftnlen diag_len; { return f06yjf_0_(0, side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb, side_len, uplo_len, transa_len, diag_len); } /* Subroutine */ int dtrsm_(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb, side_len, uplo_len, transa_len, diag_len) char *side, *uplo, *transa, *diag; integer *m, *n; doublereal *alpha, *a; integer *lda; doublereal *b; integer *ldb; ftnlen side_len; ftnlen uplo_len; ftnlen transa_len; ftnlen diag_len; { return f06yjf_0_(1, side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb, side_len, uplo_len, transa_len, diag_len); } /* Subroutine */ int f07adg_(m, n, a, lda, piv, info) integer *m, *n; doublereal *a; integer *lda; doublereal *piv; integer *info; { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5; doublereal d__1; /* Builtin functions */ double sqrt(); /* Local variables */ static integer i__, j; extern /* Subroutine */ int f07adh_(), f07adj_(), f06aaz_(), dgemm_(); static integer iinfo; extern /* Subroutine */ int f07zaz_(), dtrsm_(); static integer jb, nb; /* MARK 15 RELEASE. NAG COPYRIGHT 1991. */ /* Purpose */ /* ======= */ /* F07ADG computes an LU factorization of a general m-by-n matrix A */ /* using partial pivoting with row interchanges. */ /* The factorization has the form */ /* A = P * L * U */ /* where P is a permutation matrix, L is lower triangular (lower */ /* trapezoidal if m > n), and U is upper triangular with unit diagonal */ /* elements (upper trapezoidal if m < n). */ /* This is the Level 3 BLAS version of the Crout algorithm. */ /* Arguments */ /* ========= */ /* M (input) INTEGER */ /* The number of rows of the matrix A. M >= 0. */ /* N (input) INTEGER */ /* The number of columns of the matrix A. N >= 0. */ /* A (input/output) REAL array, dimension (LDA,N) */ /* On entry, the m by n matrix to be factored. */ /* On exit, the factors L and U from the factorization */ /* A = P*L*U; the unit diagonal elements of U are not stored. */ /* LDA (input) INTEGER */ /* The leading dimension of the array A. LDA >= max(1,M). */ /* PIV (output) REAL array, dimension (M) */ /* The pivot indices; for 1 <= i <= min(M,N), row i of the */ /* matrix was interchanged with row PIV(i). The rest of PIV is */ /* used for workspace. */ /* INFO (output) INTEGER */ /* = 0: successful exit */ /* < 0: if INFO = -k, the k-th argument had an illegal value */ /* > 0: if INFO = k, approximate singularity has been detected */ /* at the k-th stage; the factorization has not been */ /* completed. */ /* This is a modified version of the LAPACK routine F07ADF/DGETRF, in */ /* which the INTEGER array IPIV has been replaced by a REAL array PIV, */ /* row-equilibration is used in the choice of pivot, U has unit diagonal */ /* elements, and the routine exits immediately if approximate */ /* singularity is detected. */ /* ===================================================================== */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Test the input parameters. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; --piv; /* Function Body */ *info = 0; if (*m < 0) { *info = -1; } else if (*n < 0) { *info = -2; } else if (*lda < max(1,*m)) { *info = -4; } if (*info != 0) { i__1 = -(*info); f06aaz_("F07ADG ", &i__1, 13L); return 0; } /* Quick return if possible */ if (*m == 0 || *n == 0) { return 0; } /* Compute 2-norm of each row and store reciprocal in PIV. */ i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { piv[i__] = 0.; /* L20: */ } i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { /* Computing 2nd power */ d__1 = a[i__ + j * a_dim1]; piv[i__] += d__1 * d__1; /* L40: */ } /* L60: */ } i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { if (piv[i__] <= 0.) { *info = i__; return 0; } else { piv[i__] = 1. / sqrt(piv[i__]); } /* L80: */ } /* Determine the block size for this environment. */ f07zaz_(&c__1, "F07ADG", &nb, &c__0, 6L); if (nb <= 1) { /* Use unblocked code. */ f07adh_(m, n, &a[a_offset], lda, &piv[1], info); } else { /* Use blocked code. */ i__1 = min(*m,*n); i__2 = nb; for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { /* Computing MIN */ i__3 = min(*m,*n) - j + 1; jb = min(i__3,nb); /* Update diagonal and subdiagonal blocks. */ i__3 = *m - j + 1; i__4 = j - 1; dgemm_("No transpose", "No transpose", &i__3, &jb, &i__4, &c_b24, &a[j + a_dim1], lda, &a[j * a_dim1 + 1], lda, &c_b15, &a[ j + j * a_dim1], lda, 12L, 12L); /* Factorize diagonal and subdiagonal blocks and test fo r */ /* approximate singularity. */ i__3 = *m - j + 1; f07adh_(&i__3, &jb, &a[j + j * a_dim1], lda, &piv[j], &iinfo); if (iinfo > 0) { *info = iinfo + j - 1; return 0; } /* Update pivot indices and apply the interchanges to co lumns */ /* 1:J-1. */ /* Computing MIN */ i__4 = *m, i__5 = j + jb - 1; i__3 = min(i__4,i__5); for (i__ = j; i__ <= i__3; ++i__) { piv[i__] = j - 1 + piv[i__]; /* L100: */ } i__3 = j - 1; i__4 = j + jb - 1; f07adj_(&i__3, &a[a_offset], lda, &j, &i__4, &piv[1], &c__1); if (j + jb <= *n) { /* Apply the interchanges to columns J+JB:N. */ i__3 = *n - j - jb + 1; i__4 = j + jb - 1; f07adj_(&i__3, &a[(j + jb) * a_dim1 + 1], lda, &j, &i__4, & piv[1], &c__1); /* Compute block row of U. */ i__3 = *n - j - jb + 1; i__4 = j - 1; dgemm_("No transpose", "No transpose", &jb, &i__3, &i__4, & c_b24, &a[j + a_dim1], lda, &a[(j + jb) * a_dim1 + 1], lda, &c_b15, &a[j + (j + jb) * a_dim1], lda, 12L, 12L); i__3 = *n - j - jb + 1; dtrsm_("Left", "Lower", "No transpose", "Non-unit", &jb, & i__3, &c_b15, &a[j + j * a_dim1], lda, &a[j + (j + jb) * a_dim1], lda, 4L, 5L, 12L, 8L); } /* L120: */ } } return 0; /* End of F07ADG */ } /* f07adg_ */ /* Subroutine */ int f07adh_(m, n, a, lda, piv, info) integer *m, *n; doublereal *a; integer *lda; doublereal *piv; integer *info; { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; doublereal d__1; /* Local variables */ static integer i__, j; extern /* Subroutine */ int dscal_(), f06aaz_(); static doublereal x, y; extern doublereal x02ajf_(); extern /* Subroutine */ int dgemv_(), dswap_(); static integer jp; static doublereal thresh; /* MARK 15 RELEASE. NAG COPYRIGHT 1991. */ /* Purpose */ /* ======= */ /* F07ADH computes an LU factorization of a general m-by-n matrix A */ /* using partial pivoting with row interchanges. */ /* The factorization has the form */ /* A = P * L * U */ /* where P is a permutation matrix, L is lower triangular (lower */ /* trapezoidal if m > n), and U is upper triangular with unit diagonal */ /* elements (upper trapezoidal if m < n). */ /* This is the Level 2 BLAS version of the Crout algorithm. */ /* Arguments */ /* ========= */ /* M (input) INTEGER */ /* The number of rows of the matrix A. M >= 0. */ /* N (input) INTEGER */ /* The number of columns of the matrix A. N >= 0. */ /* A (input/output) REAL array, dimension (LDA,N) */ /* On entry, the m by n matrix to be factored. */ /* On exit, the factors L and U from the factorization */ /* A = P*L*U; the unit diagonal elements of U are not stored. */ /* LDA (input) INTEGER */ /* The leading dimension of the array A. LDA >= max(1,M). */ /* PIV (input/output) REAL array, dimension (M) */ /* On entry, M scale factors for equilibrating the rows of A. */ /* On exit, the pivot indices; for 1 <= i <= min(M,N), row i of */ /* the matrix was interchanged with row PIV(i). */ /* INFO (output) INTEGER */ /* = 0: successful exit */ /* < 0: if INFO = -k, the k-th argument had an illegal value */ /* > 0: if INFO = k, approximate singularity has been detected a */ /* the k-th stage; the factorization has not been */ /* completed. */ /* This is a modified version of the LAPACK routine F07ADZ/DGETF2, in */ /* which the INTEGER array IPIV has been replaced by a REAL array PIV, */ /* row-equilibration is used in the choice of pivot, U has unit diagonal */ /* elements, and the routine exits immediately if approximate */ /* singularity is detected. */ /* ===================================================================== */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Test the input parameters. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; --piv; /* Function Body */ *info = 0; if (*m < 0) { *info = -1; } else if (*n < 0) { *info = -2; } else if (*lda < max(1,*m)) { *info = -4; } if (*info != 0) { i__1 = -(*info); f06aaz_("F07ADH ", &i__1, 13L); return 0; } /* Quick return if possible */ if (*m == 0 || *n == 0) { return 0; } /* Set threshold for test for singularity */ thresh = x02ajf_() * 8.; i__1 = min(*m,*n); for (j = 1; j <= i__1; ++j) { /* Update diagonal and subdiagonal elements in column J. */ i__2 = *m - j + 1; i__3 = j - 1; dgemv_("No transpose", &i__2, &i__3, &c_b24, &a[j + a_dim1], lda, &a[ j * a_dim1 + 1], &c__1, &c_b15, &a[j + j * a_dim1], &c__1, 12L); /* Find pivot and test for singularity. */ jp = j; x = 0.; i__2 = *m; for (i__ = j; i__ <= i__2; ++i__) { y = (d__1 = a[i__ + j * a_dim1], abs(d__1)) * piv[i__]; if (y > x) { jp = i__; x = y; } /* L20: */ } piv[jp] = piv[j]; piv[j] = (doublereal) jp; if (x >= thresh) { /* Apply interchange to columns 1:N. */ if (jp != j) { dswap_(n, &a[j + a_dim1], lda, &a[jp + a_dim1], lda); } if (j + 1 <= *n) { /* Compute row of U. */ i__2 = j - 1; i__3 = *n - j; dgemv_("Transpose", &i__2, &i__3, &c_b24, &a[(j + 1) * a_dim1 + 1], lda, &a[j + a_dim1], lda, &c_b15, &a[j + (j + 1) * a_dim1], lda, 9L); i__2 = *n - j; d__1 = 1. / a[j + j * a_dim1]; dscal_(&i__2, &d__1, &a[j + (j + 1) * a_dim1], lda); } } else { /* If A( JP, J ) is small, set INFO to indicate that a s mall */ /* pivot has been found. */ *info = j; return 0; } /* L40: */ } return 0; /* End of F07ADH */ } /* f07adh_ */ /* Subroutine */ int f07adj_(n, a, lda, k1, k2, piv, incx) integer *n; doublereal *a; integer *lda, *k1, *k2; doublereal *piv; integer *incx; { /* System generated locals */ integer a_dim1, a_offset, i__1; /* Builtin functions */ integer i_dnnt(); /* Local variables */ static integer i__; extern /* Subroutine */ int dswap_(); static integer ip, ix; /* MARK 15 RELEASE. NAG COPYRIGHT 1991. */ /* Purpose */ /* ======= */ /* F07ADJ performs a series of row interchanges on the matrix A. */ /* One row interchange is initiated for each of rows K1 through K2 of A. */ /* Arguments */ /* ========= */ /* N (input) INTEGER */ /* The number of columns of the matrix A. */ /* A (input/output) REAL array, dimension (LDA,N) */ /* On entry, the matrix of column dimension N to which the row */ /* interchanges will be applied. */ /* On exit, the permuted matrix. */ /* LDA (input) INTEGER */ /* The leading dimension of the array A. */ /* K1 (input) INTEGER */ /* The first element of PIV for which a row interchange will */ /* be done. */ /* K2 (input) INTEGER */ /* The last element of PIV for which a row interchange will */ /* be done. */ /* PIV (input) REAL array, dimension( M*abs(INCX) ) */ /* The vector of pivot indices. Only the elements in positions */ /* K1 through K2 of PIV are accessed. */ /* PIV(K) = L implies rows K and L are to be interchanged. */ /* INCX (input) INTEGER */ /* The increment between succesive values of PIV. If PIV */ /* is negative, the pivots are applied in reverse order. */ /* This is a modified version of the LAPACK auxiliary routine */ /* F07ADY/DLASWP, in which the INTEGER array IPIV has been replaced by */ /* a REAL array PIV. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Interchange row I with row PIV(I) for each of rows K1 through K2. */ /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; --piv; /* Function Body */ if (*incx == 0) { return 0; } if (*incx > 0) { ix = *k1; } else { ix = (1 - *k2) * *incx + 1; } if (*incx == 1) { i__1 = *k2; for (i__ = *k1; i__ <= i__1; ++i__) { ip = i_dnnt(&piv[i__]); if (ip != i__) { dswap_(n, &a[i__ + a_dim1], lda, &a[ip + a_dim1], lda); } /* L20: */ } } else if (*incx > 1) { i__1 = *k2; for (i__ = *k1; i__ <= i__1; ++i__) { ip = i_dnnt(&piv[ix]); if (ip != i__) { dswap_(n, &a[i__ + a_dim1], lda, &a[ip + a_dim1], lda); } ix += *incx; /* L40: */ } } else if (*incx < 0) { i__1 = *k1; for (i__ = *k2; i__ >= i__1; --i__) { ip = i_dnnt(&piv[ix]); if (ip != i__) { dswap_(n, &a[i__ + a_dim1], lda, &a[ip + a_dim1], lda); } ix += *incx; /* L60: */ } } return 0; /* End of F07ADJ */ } /* f07adj_ */ integer f07zay_(name__, name_len) char *name__; ftnlen name_len; { /* System generated locals */ integer ret_val; /* Local variables */ static char name4[1], name5[1]; static integer j, k; /* MARK 15 RELEASE. NAG COPYRIGHT 1991. */ /* F07ZAY returns a unique positive integer code */ /* corresponding to a six-letter NAG routine name */ /* given in NAME. If NAME is not recognised, 0 */ /* is returned. */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Executable Statements .. */ if (*(unsigned char *)&name__[2] == '7') { *(unsigned char *)name4 = *(unsigned char *)&name__[3]; *(unsigned char *)name5 = *(unsigned char *)&name__[4]; } else { *(unsigned char *)name4 = *(unsigned char *)name__; *(unsigned char *)name5 = *(unsigned char *)&name__[1]; } if (*(unsigned char *)name4 == 'A') { j = 0; } else if (*(unsigned char *)name4 == 'B') { j = 1; } else if (*(unsigned char *)name4 == 'F') { j = 2; } else if (*(unsigned char *)name4 == 'H') { j = 3; } else if (*(unsigned char *)name4 == 'M') { j = 4; } else if (*(unsigned char *)name4 == 'N') { j = 5; } else if (*(unsigned char *)name4 == 'T') { j = 6; } else { j = -1; } if (*(unsigned char *)name5 == 'D') { k = 0; } else if (*(unsigned char *)name5 == 'J') { k = 1; } else if (*(unsigned char *)name5 == 'R') { k = 2; } else if (*(unsigned char *)name5 == 'W') { k = 3; } else { k = -1; } if (j < 0 || k < 0) { ret_val = 0; } else { /* F07ZAY is in the range 1-28. */ ret_val = (j << 2) + 1 + k; } return ret_val; } /* f07zay_ */ /* Subroutine */ int f07zaz_(ispec, name__, ival, rwflag, name_len) integer *ispec; char *name__; integer *ival, *rwflag; ftnlen name_len; { /* Initialized data */ static integer iparms[51] /* was [3][17] */ = { 1,0,0,1,1,0,1,0,0,1,1,0, 1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,80,25,0,48,10,0,1, 1,0,1,0,0,1,0,0 }; static integer point[28] = { 1,2,3,4,5,0,6,0,7,8,9,10,11,0,12,0,13,0,14,0, 0,0,15,0,0,16,0,17 }; static integer icode; extern integer f07zay_(); /* Mark 15 Release. NAG Copyright 1991 */ /* -- NAG version of LAPACK auxiliary routine ILAENV */ /* This version generated by program GENZAZ */ /* Purpose */ /* ======= */ /* F07ZAZ sets or returns problem-dependent */ /* parameters for the local environment. See */ /* ISPEC for a description of the parameters. */ /* The problem-dependent parameters are contained */ /* in the integer array IPARMS, and the value with */ /* index ISPEC is set or copied to IVAL. */ /* Arguments */ /* ========= */ /* ISPEC (input) INTEGER */ /* Specifies the parameter to be set or */ /* returned by F07ZAZ. */ /* = 1: the optimal blocksize; if this value */ /* is 1, an unblocked algorithm will give */ /* the best performance. */ /* = 2: the minimum block size for which the */ /* block routine should be used; if the */ /* usable block size is less than this */ /* value, an unblocked routine should be */ /* used. */ /* = 3: the crossover point (for N less than */ /* this value, an unblocked routine should */ /* be used) */ /* NAME (input) CHARACTER*(*) */ /* The name of the calling subroutine. */ /* IVAL (input/output) INTEGER */ /* the value of the parameter set or returned. */ /* FLAG (input) INTEGER */ /* = 0: F07ZAZ returns in IVAL the value of */ /* the parameter specified by ISPEC. */ /* = 1: F07ZAZ sets the parameter specified */ /* by ISPEC to the value in IVAL. */ /* ============================================== */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. Save statement .. */ /* .. Data statements .. */ /* .. Executable Statements .. */ /* Convert the NAG name to an integer code. */ icode = f07zay_(name__, name_len); if (*ispec < 1 || *ispec > 3) { /* Invalid value for ISPEC */ *ival = -1; } else if (icode == 0) { /* Invalid value for NAME */ *ival = -2; } else if (point[icode - 1] == 0) { /* Invalid value for NAME */ *ival = -2; } else if (*rwflag == 0) { /* Read the value of a parameter */ *ival = iparms[*ispec + point[icode - 1] * 3 - 4]; } else { /* Set the value of a parameter */ iparms[*ispec + point[icode - 1] * 3 - 4] = *ival; } return 0; /* End of F07ZAZ */ } /* f07zaz_ */ integer p01abf_(ifail, ierror, srname, nrec, rec, srname_len, rec_len) integer *ifail, *ierror; char *srname; integer *nrec; char *rec; ftnlen srname_len; ftnlen rec_len; { /* Format strings */ static char fmt_99999[] = "(\002 ** ABNORMAL EXIT from NAG Library routi\ ne \002,a,\002: IFAIL\002,\002 =\002,i6)"; /* System generated locals */ integer ret_val, i__1; /* Builtin functions */ integer s_wsfi(), do_fio(), e_wsfi(); /* Local variables */ static integer nerr; static char mess[72]; static integer i__; extern /* Subroutine */ int x04aaf_(), x04baf_(), p01abz_(); /* Fortran I/O blocks */ static icilist io___330 = { 0, mess, 0, fmt_99999, 72, 1 }; /* MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. */ /* MARK 13 REVISED. IER-621 (APR 1988). */ /* MARK 13B REVISED. IER-668 (AUG 1988). */ /* P01ABF is the error-handling routine for the NAG Library. */ /* P01ABF either returns the value of IERROR through the routine */ /* name (soft failure), or terminates execution of the program */ /* (hard failure). Diagnostic messages may be output. */ /* If IERROR = 0 (successful exit from the calling routine), */ /* the value 0 is returned through the routine name, and no */ /* message is output */ /* If IERROR is non-zero (abnormal exit from the calling routine), */ /* the action taken depends on the value of IFAIL. */ /* IFAIL = 1: soft failure, silent exit (i.e. no messages are */ /* output) */ /* IFAIL = -1: soft failure, noisy exit (i.e. messages are output) */ /* IFAIL =-13: soft failure, noisy exit but standard messages from */ /* P01ABF are suppressed */ /* IFAIL = 0: hard failure, noisy exit */ /* For compatibility with certain routines included before Mark 12 */ /* P01ABF also allows an alternative specification of IFAIL in which */ /* it is regarded as a decimal integer with least significant digits */ /* cba. Then */ /* a = 0: hard failure a = 1: soft failure */ /* b = 0: silent exit b = 1: noisy exit */ /* except that hard failure now always implies a noisy exit. */ /* S.Hammarling, M.P.Hooper and J.J.du Croz, NAG Central Office. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ rec -= rec_len; /* Function Body */ if (*ierror != 0) { /* Abnormal exit from calling routine */ if (*ifail == -1 || *ifail == 0 || *ifail == -13 || *ifail > 0 && * ifail / 10 % 10 != 0) { /* Noisy exit */ x04aaf_(&c__0, &nerr); i__1 = *nrec; for (i__ = 1; i__ <= i__1; ++i__) { x04baf_(&nerr, rec + i__ * rec_len, rec_len); /* L20: */ } if (*ifail != -13) { s_wsfi(&io___330); do_fio(&c__1, srname, srname_len); do_fio(&c__1, (char *)&(*ierror), (ftnlen)sizeof(integer)); e_wsfi(); x04baf_(&nerr, mess, 72L); if ((i__1 = *ifail % 10, abs(i__1)) != 1) { /* Hard failure */ x04baf_(&nerr, " ** NAG hard failure - execution termina\ ted", 43L); p01abz_(); } else { /* Soft failure */ x04baf_(&nerr, " ** NAG soft failure - control returned", 39L); } } } } ret_val = *ierror; return ret_val; } /* p01abf_ */ /* Subroutine */ int p01abw_(n, name__, inform__, ierr, srname, n_len, name_len, srname_len) char *n, *name__; integer *inform__, *ierr; char *srname; ftnlen n_len; ftnlen name_len; ftnlen srname_len; { /* Format strings */ static char fmt_99999[] = "(\002 ***** Parameter \002,a,\002 is inval\ id in routine \002,a,\002 ***** \002,/8x,\002Value supplied is\002,/8x,a)"; /* Builtin functions */ integer s_wsfi(), do_fio(), e_wsfi(); /* Local variables */ static integer nerr; extern /* Subroutine */ int x04aaf_(), x04baf_(); static char rec[65*3]; /* Fortran I/O blocks */ static icilist io___333 = { 0, rec, 0, fmt_99999, 65, 3 }; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* P01ABW increases the value of IERR by 1 and, if */ /* ( mod( INFORM, 10 ).ne.1 ).or.( mod( INFORM/10, 10 ).ne.0 ) */ /* writes a message on the current error message channel giving the */ /* value of N, a message to say that N is invalid and the strings */ /* NAME and SRNAME. */ /* NAME must be the name of the actual argument for N and SRNAME must */ /* be the name of the calling routine. */ /* This routine is intended for use when N is an invalid input */ /* parameter to routine SRNAME. For example */ /* IERR = 0 */ /* IF( N.NE.'Valid value' ) */ /* $ CALL P01ABW( N, 'N', IDIAG, IERR, SRNAME ) */ /* -- Written on 15-November-1984. */ /* Sven Hammarling, Nag Central Office. */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ ++(*ierr); if (*inform__ % 10 != 1 || *inform__ / 10 % 10 != 0) { x04aaf_(&c__0, &nerr); s_wsfi(&io___333); do_fio(&c__1, name__, name_len); do_fio(&c__1, srname, srname_len); do_fio(&c__1, n, n_len); e_wsfi(); x04baf_(&nerr, " ", 1L); x04baf_(&nerr, rec, 65L); x04baf_(&nerr, rec + 65, 65L); x04baf_(&nerr, rec + 130, 65L); } return 0; /* End of P01ABW. */ } /* p01abw_ */ /* Subroutine */ int p01aby_(n, name__, inform__, ierr, srname, name_len, srname_len) integer *n; char *name__; integer *inform__, *ierr; char *srname; ftnlen name_len; ftnlen srname_len; { /* Format strings */ static char fmt_99999[] = "(\002 ***** Parameter \002,a,\002 is inval\ id in routine \002,a,\002 ***** \002,/8x,\002Value supplied is \002,i6)"; /* Builtin functions */ integer s_wsfi(), do_fio(), e_wsfi(); /* Local variables */ static integer nerr; extern /* Subroutine */ int x04aaf_(), x04baf_(); static char rec[65*2]; /* Fortran I/O blocks */ static icilist io___336 = { 0, rec, 0, fmt_99999, 65, 2 }; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* P01ABY increases the value of IERR by 1 and, if */ /* ( mod( INFORM, 10 ).ne.1 ).or.( mod( INFORM/10, 10 ).ne.0 ) */ /* writes a message on the current error message channel giving the */ /* value of N, a message to say that N is invalid and the strings */ /* NAME and SRNAME. */ /* NAME must be the name of the actual argument for N and SRNAME must */ /* be the name of the calling routine. */ /* This routine is intended for use when N is an invalid input */ /* parameter to routine SRNAME. For example */ /* IERR = 0 */ /* IF( N.LT.1 )CALL P01ABY( N, 'N', IDIAG, IERR, SRNAME ) */ /* -- Written on 23-February-1984. Sven. */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ ++(*ierr); if (*inform__ % 10 != 1 || *inform__ / 10 % 10 != 0) { x04aaf_(&c__0, &nerr); s_wsfi(&io___336); do_fio(&c__1, name__, name_len); do_fio(&c__1, srname, srname_len); do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer)); e_wsfi(); x04baf_(&nerr, " ", 1L); x04baf_(&nerr, rec, 65L); x04baf_(&nerr, rec + 65, 65L); } return 0; /* End of P01ABY. */ } /* p01aby_ */ /* Subroutine */ int p01abz_() { /* Builtin functions */ /* Subroutine */ int s_stop(); /* MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. */ /* Terminates execution when a hard failure occurs. */ /* ******************** IMPLEMENTATION NOTE ******************** */ /* The following STOP statement may be replaced by a call to an */ /* implementation-dependent routine to display a message and/or */ /* to abort the program. */ /* ************************************************************* */ /* .. Executable Statements .. */ /* s_stop("", 0L); */ return(0); } /* p01abz_ */ integer p01acf_(ifail, ierror, srname, varbnm, nrec, rec, srname_len, varbnm_len, rec_len) integer *ifail, *ierror; char *srname, *varbnm; integer *nrec; char *rec; ftnlen srname_len; ftnlen varbnm_len; ftnlen rec_len; { /* Format strings */ static char fmt_99999[] = "(\002 ** ABNORMAL EXIT from NAG Library routi\ ne \002,a,\002: \002,a,\002 =\002,i6)"; static char fmt_99998[] = "(\002 ** ABNORMAL EXIT from NAG Library routi\ ne \002,a)"; /* System generated locals */ integer ret_val, i__1; /* Builtin functions */ integer i_len(), s_wsfi(), do_fio(), e_wsfi(); /* Local variables */ static integer nerr; static char mess[72]; static integer i__; extern /* Subroutine */ int x04aaf_(), x04baf_(), p01abz_(); static integer varlen; /* Fortran I/O blocks */ static icilist io___341 = { 0, mess, 0, fmt_99999, 72, 1 }; static icilist io___342 = { 0, mess, 0, fmt_99998, 72, 1 }; /* MARK 15 RELEASE. NAG COPYRIGHT 1991. */ /* P01ACF is the error-handling routine for the F06 AND F07 */ /* Chapters of the NAG Fortran Library. It is a slightly modified */ /* version of P01ABF. */ /* P01ACF either returns the value of IERROR through the routine */ /* name (soft failure), or terminates execution of the program */ /* (hard failure). Diagnostic messages may be output. */ /* If IERROR = 0 (successful exit from the calling routine), */ /* the value 0 is returned through the routine name, and no */ /* message is output */ /* If IERROR is non-zero (abnormal exit from the calling routine), */ /* the action taken depends on the value of IFAIL. */ /* IFAIL = 1: soft failure, silent exit (i.e. no messages are */ /* output) */ /* IFAIL = -1: soft failure, noisy exit (i.e. messages are output) */ /* IFAIL =-13: soft failure, noisy exit but standard messages from */ /* P01ACF are suppressed */ /* IFAIL = 0: hard failure, noisy exit */ /* For compatibility with certain routines included before Mark 12 */ /* P01ACF also allows an alternative specification of IFAIL in which */ /* it is regarded as a decimal integer with least significant digits */ /* cba. Then */ /* a = 0: hard failure a = 1: soft failure */ /* b = 0: silent exit b = 1: noisy exit */ /* except that hard failure now always implies a noisy exit. */ /* S.Hammarling, M.P.Hooper and J.J.du Croz, NAG Central Office. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. External Subroutines .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ rec -= rec_len; /* Function Body */ if (*ierror != 0) { varlen = 0; for (i__ = i_len(varbnm, varbnm_len); i__ >= 1; --i__) { if (*(unsigned char *)&varbnm[i__ - 1] != ' ') { varlen = i__; goto L40; } /* L20: */ } L40: /* Abnormal exit from calling routine */ if (*ifail == -1 || *ifail == 0 || *ifail == -13 || *ifail > 0 && * ifail / 10 % 10 != 0) { /* Noisy exit */ x04aaf_(&c__0, &nerr); i__1 = *nrec; for (i__ = 1; i__ <= i__1; ++i__) { x04baf_(&nerr, rec + i__ * rec_len, rec_len); /* L60: */ } if (*ifail != -13) { if (varlen != 0) { s_wsfi(&io___341); do_fio(&c__1, srname, srname_len); do_fio(&c__1, varbnm, varlen); do_fio(&c__1, (char *)&(*ierror), (ftnlen)sizeof(integer)) ; e_wsfi(); } else { s_wsfi(&io___342); do_fio(&c__1, srname, srname_len); e_wsfi(); } x04baf_(&nerr, mess, 72L); if ((i__1 = *ifail % 10, abs(i__1)) != 1) { /* Hard failure */ x04baf_(&nerr, " ** NAG hard failure - execution termina\ ted", 43L); p01abz_(); } else { /* Soft failure */ x04baf_(&nerr, " ** NAG soft failure - control returned", 39L); } } } } ret_val = *ierror; return ret_val; } /* p01acf_ */ doublereal x02ajf_() { /* Initialized data */ static doublereal x02con = 1.11022302462516e-16; /* System generated locals */ doublereal ret_val; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* RETURNS (1/2)*B**(1-P) IF ROUNDS IS .TRUE. */ /* RETURNS B**(1-P) OTHERWISE */ /* .. Executable Statements .. */ ret_val = x02con; return ret_val; } /* x02ajf_ */ doublereal x02akf_() { /* Initialized data */ static doublereal x02con = 2.22507385850721e-308; /* System generated locals */ doublereal ret_val; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* RETURNS B**(EMIN-1) (THE SMALLEST POSITIVE MODEL NUMBER) */ /* .. Executable Statements .. */ ret_val = x02con; return ret_val; } /* x02akf_ */ doublereal x02amf_() { /* Initialized data */ static doublereal x02con = 2.22507385850739e-308; /* System generated locals */ doublereal ret_val; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* RETURNS THE 'SAFE RANGE' PARAMETER */ /* I.E. THE SMALLEST POSITIVE MODEL NUMBER Z SUCH THAT */ /* FOR ANY X WHICH SATISFIES X.GE.Z AND X.LE.1/Z */ /* THE FOLLOWING CAN BE COMPUTED WITHOUT OVERFLOW, UNDERFLOW OR OTHER */ /* ERROR */ /* -X */ /* 1.0/X */ /* SQRT(X) */ /* LOG(X) */ /* EXP(LOG(X)) */ /* Y**(LOG(X)/LOG(Y)) FOR ANY Y */ /* .. Executable Statements .. */ ret_val = x02con; return ret_val; } /* x02amf_ */ integer x02bhf_() { /* System generated locals */ integer ret_val; /* MARK 12 RELEASE. NAG COPYRIGHT 1986. */ /* RETURNS THE MODEL PARAMETER, B. */ /* .. Executable Statements .. */ ret_val = 2; return ret_val; } /* x02bhf_ */ /* Subroutine */ int x03aaf_(a, isizea, b, isizeb, n, istepa, istepb, c1, c2, d1, d2, sw, ifail) doublereal *a; integer *isizea; doublereal *b; integer *isizeb, *n, *istepa, *istepb; doublereal *c1, *c2, *d1, *d2; logical *sw; integer *ifail; { /* System generated locals */ integer i__1; /* Local variables */ static integer ierr, i__; extern integer p01abf_(); static char p01rec[1*1]; extern /* Subroutine */ int x03aay_(); static integer is, it; static doublereal sum; /* MARK 10 RE-ISSUE. NAG COPYRIGHT 1982 */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 12 REVISED. IER-524 (AUG 1986). */ /* DOUBLE PRECISION BASE VERSION */ /* CALCULATES THE VALUE OF A SCALAR PRODUCT USING BASIC OR */ /* ADDITIONAL PRECISION AND ADDS IT TO A BASIC OR ADDITIONAL */ /* PRECISION INITIAL VALUE. */ /* FOR THIS DOUBLE PRECISION VERSION, ALL ADDITIONAL (I.E. */ /* QUADRUPLE) PRECISION COMPUTATION IS PERFORMED BY THE AUXILIARY */ /* ROUTINE X03AAY. SEE THE COMMENTS AT THE HEAD OF THAT ROUTINE */ /* CONCERNING IMPLEMENTATION. */ /* .. Parameters .. */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Local Arrays .. */ /* .. External Functions .. */ /* .. External Subroutines .. */ /* .. Executable Statements .. */ /* Parameter adjustments */ --a; --b; /* Function Body */ ierr = 0; if (*istepa <= 0 || *istepb <= 0) { ierr = 1; } if (*isizea < (*n - 1) * *istepa + 1 || *isizeb < (*n - 1) * *istepb + 1) { ierr = 2; } if (ierr == 0) { goto L20; } *ifail = p01abf_(ifail, &ierr, "X03AAF", &c__0, p01rec, 6L, 1L); return 0; L20: *ifail = 0; if (*sw) { goto L80; } /* BASIC PRECISION CALCULATION */ sum = *c1; if (*n < 1) { goto L60; } is = 1; it = 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { sum += a[is] * b[it]; is += *istepa; it += *istepb; /* L40: */ } L60: *d1 = sum; *d2 = 0.; return 0; /* ADDITIONAL PRECISION COMPUTATION */ L80: x03aay_(&a[1], isizea, &b[1], isizeb, n, istepa, istepb, c1, c2, d1, d2); return 0; } /* x03aaf_ */ /* Subroutine */ int x03aay_(a, isizea, b, isizeb, n, istepa, istepb, c1, c2, d1, d2) doublereal *a; integer *isizea; doublereal *b; integer *isizeb, *n, *istepa, *istepb; doublereal *c1, *c2, *d1, *d2; { /* Initialized data */ static doublereal cons = 134217729.; /* System generated locals */ integer i__1; /* Local variables */ static doublereal summ; static integer i__; static doublereal r__, s, z__, aa, bb, ha, hb, ta, tb; static integer is, it; static doublereal zz, sum; /* MARK 10 RELEASE. NAG COPYRIGHT 1982 */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* DOUBLE PRECISION BASE VERSION */ /* PERFORMS QUADRUPLE PRECISION COMPUTATION FOR X03AAF */ /* ************************************************************** */ /* THIS FORTRAN CODE WILL NOT WORK ON ALL MACHINES. IT IS */ /* PRESUMED TO WORK IF THE MACHINE SATISFIES ONE OF THE FOLLOWING */ /* ASSUMPTIONS. */ /* A.) THERE ARE AN EVEN NUMBER OF B-ARY DIGITS IN THE MANTISSA */ /* - OF A DOUBLE PRECISION NUMBER (WHERE B IS THE BASE FOR THE */ /* - REPRESENTATION OF FLOATING-POINT NUMBERS), AND THE */ /* - COMPUTED RESULT OF A DOUBLE PRECISION ADDITION, */ /* - SUBTRACTION OR MULTIPLICATION IS EITHER CORRECTLY ROUNDED */ /* - OR CORRECTLY CHOPPED. */ /* B.) FLOATING-POINT NUMBERS ARE REPRESENTED TO THE BASE 2 (WITH */ /* - ANY NUMBER OF BITS IN THE MANTISSA OF A DOUBLE PRECISION */ /* - NUMBER), AND THE COMPUTED RESULT OF A DOUBLE PRECISION */ /* - ADDITION, SUBTRACTION OR MULTIPLICATION IS CORRECTLY */ /* - ROUNDED. */ /* REFERENCES- */ /* T.J. DEKKER A FLOATING-POINT TECHNIQUE FOR EXTENDING THE */ /* AVAILABLE PRECISION. NUMER. MATH. 18, 224-242 (1971) */ /* S. LINNAINMAA SOFTWARE FOR DOUBLED-PRECISION FLOATING-POINT */ /* COMPUTATIONS. ACM TRANS. MATH. SOFTWARE 7, 272-283 (1981) */ /* IF THE ABOVE ASSUMPTIONS ARE NOT SATISFIED, THIS ROUTINE MUST */ /* BE IMPLEMENTED IN ASSEMBLY LANGUAGE. IN ANY CASE ASSEMBLY */ /* LANGUAGE MAY BE PREFERABLE FOR GREATER EFFICIENCY. CONSULT */ /* NAG CENTRAL OFFICE. */ /* THE ROUTINE MUST SIMULATE THE FOLLOWING QUADRUPLE PRECISION */ /* CODING IN PSEUDO-FORTRAN, WHERE */ /* - QEXTD CONVERTS FROM DOUBLE TO QUADRUPLE PRECISION */ /* - DBLEQ CONVERTS FROM QUADRUPLE TO DOUBLE PRECISION. */ /* QUADRUPLE PRECISION SUM */ /* SUM = QEXTD(C1) + QEXTD(C2) */ /* IF (N.LT.1) GO TO 80 */ /* IS = 1 */ /* IT = 1 */ /* DO 60 I = 1, N */ /* SUM = SUM + QEXTD(A(IS))*QEXTD(B(IT)) */ /* IS = IS + ISTEPA */ /* IT = IT + ISTEPB */ /* 60 CONTINUE */ /* 80 D1 = SUM -- ROUNDED TO DOUBLE PRECISION */ /* D2 = DBLEQ(SUM-QEXTD(D1)) */ /* ************************************************************** */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. Local Scalars .. */ /* .. Intrinsic Functions .. */ /* .. Data statements .. */ /* ************* IMPLEMENTATION-DEPENDENT CONSTANT ************** */ /* CONS MUST BE SET TO B**(T - INT(T/2)) + 1 , WHERE T IS THE */ /* NUMBER OF B-ARY DIGITS IN THE MANTISSA OF A DOUBLE PRECISION */ /* NUMBER. */ /* FOR B = 16 AND T = 14 (E.G. IBM 370) OR */ /* FOR B = 2 AND T = 56 (E.G. DEC VAX-11) */ /* DATA CONS /268435457.0D0/ */ /* Parameter adjustments */ --a; --b; /* Function Body */ /* ************************************************************** */ /* .. Executable Statements .. */ sum = *c1 + *c2; summ = *c1 - sum + *c2; if (*n < 1) { goto L80; } is = 1; it = 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { aa = a[is]; bb = b[it]; z__ = aa * cons; ha = aa - z__ + z__; ta = aa - ha; z__ = bb * cons; hb = bb - z__ + z__; tb = bb - hb; z__ = aa * bb; zz = ha * hb - z__ + ha * tb + ta * hb + ta * tb; r__ = z__ + sum; if (abs(z__) > abs(sum)) { goto L20; } s = sum - r__ + z__ + zz + summ; goto L40; L20: s = z__ - r__ + sum + summ + zz; L40: sum = r__ + s; summ = r__ - sum + s; is += *istepa; it += *istepb; /* L60: */ } /* 80 D1 = SUM + (SUMM+SUMM) */ /* *************** IMPLEMENTATION DEPENDENT CODE **************** */ /* THE PREVIOUS STATEMENT ASSUMES THAT THE RESULT OF A DOUBLE */ /* PRECISION ADDITION IS TRUNCATED. IF IT IS ROUNDED, THEN */ /* THE STATEMENT MUST BE CHANGED TO */ L80: *d1 = sum + summ; /* ************************************************************** */ *d2 = sum - *d1 + summ; return 0; } /* x03aay_ */ /* Subroutine */ int x04aaf_(i__, nerr) integer *i__, *nerr; { /* Initialized data */ static integer nerr1 = 6; /* MARK 7 RELEASE. NAG COPYRIGHT 1978 */ /* MARK 7C REVISED IER-190 (MAY 1979) */ /* MARK 11.5(F77) REVISED. (SEPT 1985.) */ /* MARK 14 REVISED. IER-829 (DEC 1989). */ /* IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER */ /* (STORED IN NERR1). */ /* IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO */ /* VALUE SPECIFIED BY NERR. */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Save statement .. */ /* .. Data statements .. */ /* .. Executable Statements .. */ if (*i__ == 0) { *nerr = nerr1; } if (*i__ == 1) { nerr1 = *nerr; } return 0; } /* x04aaf_ */ /* Subroutine */ int x04baf_(nout, rec, rec_len) integer *nout; char *rec; ftnlen rec_len; { /* Format strings */ static char fmt_99999[] = "(a)"; /* Builtin functions */ integer i_len(), s_wsfe(), do_fio(), e_wsfe(); /* Local variables */ static integer i__; /* Fortran I/O blocks */ static cilist io___370 = { 0, 0, 0, fmt_99999, 0 }; /* MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. */ /* X04BAF writes the contents of REC to the unit defined by NOUT. */ /* Trailing blanks are not output, except that if REC is entirely */ /* blank, a single blank character is output. */ /* If NOUT.lt.0, i.e. if NOUT is not a valid Fortran unit identifier, */ /* then no output occurs. */ /* .. Scalar Arguments .. */ /* .. Local Scalars .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ if (*nout >= 0) { /* Remove trailing blanks */ for (i__ = i_len(rec, rec_len); i__ >= 2; --i__) { if (*(unsigned char *)&rec[i__ - 1] != ' ') { goto L40; } /* L20: */ } /* Write record to external file */ L40: io___370.ciunit = *nout; s_wsfe(&io___370); do_fio(&c__1, rec, i__); e_wsfe(); } return 0; } /* x04baf_ */