From 62af5826eac61e30146b6df13fb56868c8c1517f Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 7 Apr 2011 15:58:01 +0000 Subject: [PATCH] Added a lapack routine --- ext/f2c_lapack/Makefile.in | 1 + ext/f2c_lapack/dlantr.c | 401 +++++++++++++++++++++++++++++++++++++ ext/lapack/Makefile.in | 1 + ext/lapack/dlantr.f | 277 +++++++++++++++++++++++++ 4 files changed, 680 insertions(+) create mode 100644 ext/f2c_lapack/dlantr.c create mode 100644 ext/lapack/dlantr.f diff --git a/ext/f2c_lapack/Makefile.in b/ext/f2c_lapack/Makefile.in index 68dff0ce8..16b3add23 100755 --- a/ext/f2c_lapack/Makefile.in +++ b/ext/f2c_lapack/Makefile.in @@ -67,6 +67,7 @@ dlabrd.o \ dlacpy.o \ dlamch.o \ dlange.o \ +dlangtr.o \ dlapy2.o \ dlarf.o \ dlarfb.o \ diff --git a/ext/f2c_lapack/dlantr.c b/ext/f2c_lapack/dlantr.c new file mode 100644 index 000000000..17a0e0724 --- /dev/null +++ b/ext/f2c_lapack/dlantr.c @@ -0,0 +1,401 @@ +/* dlantr.f -- translated by f2c (version 20031025). + You must link the resulting object file with libf2c: + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip +*/ + +#include "f2c.h" + +/* Table of constant values */ + +static integer c__1 = 1; + +doublereal dlantr_(char *norm, char *uplo, char *diag, integer *m, integer *n, + doublereal *a, integer *lda, doublereal *work, ftnlen norm_len, + ftnlen uplo_len, ftnlen diag_len) +{ + /* System generated locals */ + integer a_dim1, a_offset, i__1, i__2, i__3, i__4; + doublereal ret_val, d__1, d__2, d__3; + + /* Builtin functions */ + double sqrt(doublereal); + + /* Local variables */ + static integer i__, j; + static doublereal sum, scale; + static logical udiag; + extern logical lsame_(char *, char *, ftnlen, ftnlen); + static doublereal value; + extern /* Subroutine */ int dlassq_(integer *, doublereal *, integer *, + doublereal *, doublereal *); + + +/* -- LAPACK auxiliary routine (version 3.0) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */ +/* Courant Institute, Argonne National Lab, and Rice University */ +/* October 31, 1992 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* Purpose */ +/* ======= */ + +/* DLANTR returns the value of the one norm, or the Frobenius norm, or */ +/* the infinity norm, or the element of largest absolute value of a */ +/* trapezoidal or triangular matrix A. */ + +/* Description */ +/* =========== */ + +/* DLANTR returns the value */ + +/* DLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm' */ +/* ( */ +/* ( norm1(A), NORM = '1', 'O' or 'o' */ +/* ( */ +/* ( normI(A), NORM = 'I' or 'i' */ +/* ( */ +/* ( normF(A), NORM = 'F', 'f', 'E' or 'e' */ + +/* where norm1 denotes the one norm of a matrix (maximum column sum), */ +/* normI denotes the infinity norm of a matrix (maximum row sum) and */ +/* normF denotes the Frobenius norm of a matrix (square root of sum of */ +/* squares). Note that max(abs(A(i,j))) is not a matrix norm. */ + +/* Arguments */ +/* ========= */ + +/* NORM (input) CHARACTER*1 */ +/* Specifies the value to be returned in DLANTR as described */ +/* above. */ + +/* UPLO (input) CHARACTER*1 */ +/* Specifies whether the matrix A is upper or lower trapezoidal. */ +/* = 'U': Upper trapezoidal */ +/* = 'L': Lower trapezoidal */ +/* Note that A is triangular instead of trapezoidal if M = N. */ + +/* DIAG (input) CHARACTER*1 */ +/* Specifies whether or not the matrix A has unit diagonal. */ +/* = 'N': Non-unit diagonal */ +/* = 'U': Unit diagonal */ + +/* M (input) INTEGER */ +/* The number of rows of the matrix A. M >= 0, and if */ +/* UPLO = 'U', M <= N. When M = 0, DLANTR is set to zero. */ + +/* N (input) INTEGER */ +/* The number of columns of the matrix A. N >= 0, and if */ +/* UPLO = 'L', N <= M. When N = 0, DLANTR is set to zero. */ + +/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */ +/* The trapezoidal matrix A (A is triangular if M = N). */ +/* If UPLO = 'U', the leading m by n upper trapezoidal part of */ +/* the array A contains the upper trapezoidal matrix, and the */ +/* strictly lower triangular part of A is not referenced. */ +/* If UPLO = 'L', the leading m by n lower trapezoidal part of */ +/* the array A contains the lower trapezoidal matrix, and the */ +/* strictly upper triangular part of A is not referenced. Note */ +/* that when DIAG = 'U', the diagonal elements of A are not */ +/* referenced and are assumed to be one. */ + +/* LDA (input) INTEGER */ +/* The leading dimension of the array A. LDA >= max(M,1). */ + +/* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), */ +/* where LWORK >= M when NORM = 'I'; otherwise, WORK is not */ +/* referenced. */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. External Subroutines .. */ +/* .. */ +/* .. External Functions .. */ +/* .. */ +/* .. Intrinsic Functions .. */ +/* .. */ +/* .. Executable Statements .. */ + + /* Parameter adjustments */ + a_dim1 = *lda; + a_offset = 1 + a_dim1; + a -= a_offset; + --work; + + /* Function Body */ + if (min(*m,*n) == 0) { + value = 0.; + } else if (lsame_(norm, "M", (ftnlen)1, (ftnlen)1)) { + +/* Find max(abs(A(i,j))). */ + + if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) { + value = 1.; + if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { + i__1 = *n; + for (j = 1; j <= i__1; ++j) { +/* Computing MIN */ + i__3 = *m, i__4 = j - 1; + i__2 = min(i__3,i__4); + for (i__ = 1; i__ <= i__2; ++i__) { +/* Computing MAX */ + d__2 = value, d__3 = (d__1 = a[i__ + j * a_dim1], abs( + d__1)); + value = max(d__2,d__3); +/* L10: */ + } +/* L20: */ + } + } else { + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = *m; + for (i__ = j + 1; i__ <= i__2; ++i__) { +/* Computing MAX */ + d__2 = value, d__3 = (d__1 = a[i__ + j * a_dim1], abs( + d__1)); + value = max(d__2,d__3); +/* L30: */ + } +/* L40: */ + } + } + } else { + value = 0.; + if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = min(*m,j); + for (i__ = 1; i__ <= i__2; ++i__) { +/* Computing MAX */ + d__2 = value, d__3 = (d__1 = a[i__ + j * a_dim1], abs( + d__1)); + value = max(d__2,d__3); +/* L50: */ + } +/* L60: */ + } + } else { + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = *m; + for (i__ = j; i__ <= i__2; ++i__) { +/* Computing MAX */ + d__2 = value, d__3 = (d__1 = a[i__ + j * a_dim1], abs( + d__1)); + value = max(d__2,d__3); +/* L70: */ + } +/* L80: */ + } + } + } + } else if (lsame_(norm, "O", (ftnlen)1, (ftnlen)1) || *(unsigned char *) + norm == '1') { + +/* Find norm1(A). */ + + value = 0.; + udiag = lsame_(diag, "U", (ftnlen)1, (ftnlen)1); + if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + if (udiag && j <= *m) { + sum = 1.; + i__2 = j - 1; + for (i__ = 1; i__ <= i__2; ++i__) { + sum += (d__1 = a[i__ + j * a_dim1], abs(d__1)); +/* L90: */ + } + } else { + sum = 0.; + i__2 = min(*m,j); + for (i__ = 1; i__ <= i__2; ++i__) { + sum += (d__1 = a[i__ + j * a_dim1], abs(d__1)); +/* L100: */ + } + } + value = max(value,sum); +/* L110: */ + } + } else { + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + if (udiag) { + sum = 1.; + i__2 = *m; + for (i__ = j + 1; i__ <= i__2; ++i__) { + sum += (d__1 = a[i__ + j * a_dim1], abs(d__1)); +/* L120: */ + } + } else { + sum = 0.; + i__2 = *m; + for (i__ = j; i__ <= i__2; ++i__) { + sum += (d__1 = a[i__ + j * a_dim1], abs(d__1)); +/* L130: */ + } + } + value = max(value,sum); +/* L140: */ + } + } + } else if (lsame_(norm, "I", (ftnlen)1, (ftnlen)1)) { + +/* Find normI(A). */ + + if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { + if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) { + i__1 = *m; + for (i__ = 1; i__ <= i__1; ++i__) { + work[i__] = 1.; +/* L150: */ + } + i__1 = *n; + for (j = 1; j <= i__1; ++j) { +/* Computing MIN */ + i__3 = *m, i__4 = j - 1; + i__2 = min(i__3,i__4); + for (i__ = 1; i__ <= i__2; ++i__) { + work[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1)); +/* L160: */ + } +/* L170: */ + } + } else { + i__1 = *m; + for (i__ = 1; i__ <= i__1; ++i__) { + work[i__] = 0.; +/* L180: */ + } + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = min(*m,j); + for (i__ = 1; i__ <= i__2; ++i__) { + work[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1)); +/* L190: */ + } +/* L200: */ + } + } + } else { + if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) { + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + work[i__] = 1.; +/* L210: */ + } + i__1 = *m; + for (i__ = *n + 1; i__ <= i__1; ++i__) { + work[i__] = 0.; +/* L220: */ + } + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = *m; + for (i__ = j + 1; i__ <= i__2; ++i__) { + work[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1)); +/* L230: */ + } +/* L240: */ + } + } else { + i__1 = *m; + for (i__ = 1; i__ <= i__1; ++i__) { + work[i__] = 0.; +/* L250: */ + } + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = *m; + for (i__ = j; i__ <= i__2; ++i__) { + work[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1)); +/* L260: */ + } +/* L270: */ + } + } + } + value = 0.; + i__1 = *m; + for (i__ = 1; i__ <= i__1; ++i__) { +/* Computing MAX */ + d__1 = value, d__2 = work[i__]; + value = max(d__1,d__2); +/* L280: */ + } + } else if (lsame_(norm, "F", (ftnlen)1, (ftnlen)1) || lsame_(norm, "E", ( + ftnlen)1, (ftnlen)1)) { + +/* Find normF(A). */ + + if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) { + if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) { + scale = 1.; + sum = (doublereal) min(*m,*n); + i__1 = *n; + for (j = 2; j <= i__1; ++j) { +/* Computing MIN */ + i__3 = *m, i__4 = j - 1; + i__2 = min(i__3,i__4); + dlassq_(&i__2, &a[j * a_dim1 + 1], &c__1, &scale, &sum); +/* L290: */ + } + } else { + scale = 0.; + sum = 1.; + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = min(*m,j); + dlassq_(&i__2, &a[j * a_dim1 + 1], &c__1, &scale, &sum); +/* L300: */ + } + } + } else { + if (lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) { + scale = 1.; + sum = (doublereal) min(*m,*n); + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = *m - j; +/* Computing MIN */ + i__3 = *m, i__4 = j + 1; + dlassq_(&i__2, &a[min(i__3,i__4) + j * a_dim1], &c__1, & + scale, &sum); +/* L310: */ + } + } else { + scale = 0.; + sum = 1.; + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = *m - j + 1; + dlassq_(&i__2, &a[j + j * a_dim1], &c__1, &scale, &sum); +/* L320: */ + } + } + } + value = scale * sqrt(sum); + } + + ret_val = value; + return ret_val; + +/* End of DLANTR */ + +} /* dlantr_ */ + diff --git a/ext/lapack/Makefile.in b/ext/lapack/Makefile.in index 50a478cca..4dd38ef02 100755 --- a/ext/lapack/Makefile.in +++ b/ext/lapack/Makefile.in @@ -38,6 +38,7 @@ dlabrd.o \ dlacpy.o \ dlamch.o \ dlange.o \ +dlangtr.o \ dlapy2.o \ dlarf.o \ dlarfb.o \ diff --git a/ext/lapack/dlantr.f b/ext/lapack/dlantr.f new file mode 100644 index 000000000..19e9b5d92 --- /dev/null +++ b/ext/lapack/dlantr.f @@ -0,0 +1,277 @@ + DOUBLE PRECISION FUNCTION DLANTR( NORM, UPLO, DIAG, M, N, A, LDA, + $ WORK ) +* +* -- LAPACK auxiliary routine (version 3.0) -- +* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +* Courant Institute, Argonne National Lab, and Rice University +* October 31, 1992 +* +* .. Scalar Arguments .. + CHARACTER DIAG, NORM, UPLO + INTEGER LDA, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DLANTR returns the value of the one norm, or the Frobenius norm, or +* the infinity norm, or the element of largest absolute value of a +* trapezoidal or triangular matrix A. +* +* Description +* =========== +* +* DLANTR returns the value +* +* DLANTR = ( max(abs(A(i,j))), NORM = 'M' or 'm' +* ( +* ( norm1(A), NORM = '1', 'O' or 'o' +* ( +* ( normI(A), NORM = 'I' or 'i' +* ( +* ( normF(A), NORM = 'F', 'f', 'E' or 'e' +* +* where norm1 denotes the one norm of a matrix (maximum column sum), +* normI denotes the infinity norm of a matrix (maximum row sum) and +* normF denotes the Frobenius norm of a matrix (square root of sum of +* squares). Note that max(abs(A(i,j))) is not a matrix norm. +* +* Arguments +* ========= +* +* NORM (input) CHARACTER*1 +* Specifies the value to be returned in DLANTR as described +* above. +* +* UPLO (input) CHARACTER*1 +* Specifies whether the matrix A is upper or lower trapezoidal. +* = 'U': Upper trapezoidal +* = 'L': Lower trapezoidal +* Note that A is triangular instead of trapezoidal if M = N. +* +* DIAG (input) CHARACTER*1 +* Specifies whether or not the matrix A has unit diagonal. +* = 'N': Non-unit diagonal +* = 'U': Unit diagonal +* +* M (input) INTEGER +* The number of rows of the matrix A. M >= 0, and if +* UPLO = 'U', M <= N. When M = 0, DLANTR is set to zero. +* +* N (input) INTEGER +* The number of columns of the matrix A. N >= 0, and if +* UPLO = 'L', N <= M. When N = 0, DLANTR is set to zero. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,N) +* The trapezoidal matrix A (A is triangular if M = N). +* If UPLO = 'U', the leading m by n upper trapezoidal part of +* the array A contains the upper trapezoidal matrix, and the +* strictly lower triangular part of A is not referenced. +* If UPLO = 'L', the leading m by n lower trapezoidal part of +* the array A contains the lower trapezoidal matrix, and the +* strictly upper triangular part of A is not referenced. Note +* that when DIAG = 'U', the diagonal elements of A are not +* referenced and are assumed to be one. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(M,1). +* +* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), +* where LWORK >= M when NORM = 'I'; otherwise, WORK is not +* referenced. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UDIAG + INTEGER I, J + DOUBLE PRECISION SCALE, SUM, VALUE +* .. +* .. External Subroutines .. + EXTERNAL DLASSQ +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN, SQRT +* .. +* .. Executable Statements .. +* + IF( MIN( M, N ).EQ.0 ) THEN + VALUE = ZERO + ELSE IF( LSAME( NORM, 'M' ) ) THEN +* +* Find max(abs(A(i,j))). +* + IF( LSAME( DIAG, 'U' ) ) THEN + VALUE = ONE + IF( LSAME( UPLO, 'U' ) ) THEN + DO 20 J = 1, N + DO 10 I = 1, MIN( M, J-1 ) + VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + 10 CONTINUE + 20 CONTINUE + ELSE + DO 40 J = 1, N + DO 30 I = J + 1, M + VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + 30 CONTINUE + 40 CONTINUE + END IF + ELSE + VALUE = ZERO + IF( LSAME( UPLO, 'U' ) ) THEN + DO 60 J = 1, N + DO 50 I = 1, MIN( M, J ) + VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + 50 CONTINUE + 60 CONTINUE + ELSE + DO 80 J = 1, N + DO 70 I = J, M + VALUE = MAX( VALUE, ABS( A( I, J ) ) ) + 70 CONTINUE + 80 CONTINUE + END IF + END IF + ELSE IF( ( LSAME( NORM, 'O' ) ) .OR. ( NORM.EQ.'1' ) ) THEN +* +* Find norm1(A). +* + VALUE = ZERO + UDIAG = LSAME( DIAG, 'U' ) + IF( LSAME( UPLO, 'U' ) ) THEN + DO 110 J = 1, N + IF( ( UDIAG ) .AND. ( J.LE.M ) ) THEN + SUM = ONE + DO 90 I = 1, J - 1 + SUM = SUM + ABS( A( I, J ) ) + 90 CONTINUE + ELSE + SUM = ZERO + DO 100 I = 1, MIN( M, J ) + SUM = SUM + ABS( A( I, J ) ) + 100 CONTINUE + END IF + VALUE = MAX( VALUE, SUM ) + 110 CONTINUE + ELSE + DO 140 J = 1, N + IF( UDIAG ) THEN + SUM = ONE + DO 120 I = J + 1, M + SUM = SUM + ABS( A( I, J ) ) + 120 CONTINUE + ELSE + SUM = ZERO + DO 130 I = J, M + SUM = SUM + ABS( A( I, J ) ) + 130 CONTINUE + END IF + VALUE = MAX( VALUE, SUM ) + 140 CONTINUE + END IF + ELSE IF( LSAME( NORM, 'I' ) ) THEN +* +* Find normI(A). +* + IF( LSAME( UPLO, 'U' ) ) THEN + IF( LSAME( DIAG, 'U' ) ) THEN + DO 150 I = 1, M + WORK( I ) = ONE + 150 CONTINUE + DO 170 J = 1, N + DO 160 I = 1, MIN( M, J-1 ) + WORK( I ) = WORK( I ) + ABS( A( I, J ) ) + 160 CONTINUE + 170 CONTINUE + ELSE + DO 180 I = 1, M + WORK( I ) = ZERO + 180 CONTINUE + DO 200 J = 1, N + DO 190 I = 1, MIN( M, J ) + WORK( I ) = WORK( I ) + ABS( A( I, J ) ) + 190 CONTINUE + 200 CONTINUE + END IF + ELSE + IF( LSAME( DIAG, 'U' ) ) THEN + DO 210 I = 1, N + WORK( I ) = ONE + 210 CONTINUE + DO 220 I = N + 1, M + WORK( I ) = ZERO + 220 CONTINUE + DO 240 J = 1, N + DO 230 I = J + 1, M + WORK( I ) = WORK( I ) + ABS( A( I, J ) ) + 230 CONTINUE + 240 CONTINUE + ELSE + DO 250 I = 1, M + WORK( I ) = ZERO + 250 CONTINUE + DO 270 J = 1, N + DO 260 I = J, M + WORK( I ) = WORK( I ) + ABS( A( I, J ) ) + 260 CONTINUE + 270 CONTINUE + END IF + END IF + VALUE = ZERO + DO 280 I = 1, M + VALUE = MAX( VALUE, WORK( I ) ) + 280 CONTINUE + ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN +* +* Find normF(A). +* + IF( LSAME( UPLO, 'U' ) ) THEN + IF( LSAME( DIAG, 'U' ) ) THEN + SCALE = ONE + SUM = MIN( M, N ) + DO 290 J = 2, N + CALL DLASSQ( MIN( M, J-1 ), A( 1, J ), 1, SCALE, SUM ) + 290 CONTINUE + ELSE + SCALE = ZERO + SUM = ONE + DO 300 J = 1, N + CALL DLASSQ( MIN( M, J ), A( 1, J ), 1, SCALE, SUM ) + 300 CONTINUE + END IF + ELSE + IF( LSAME( DIAG, 'U' ) ) THEN + SCALE = ONE + SUM = MIN( M, N ) + DO 310 J = 1, N + CALL DLASSQ( M-J, A( MIN( M, J+1 ), J ), 1, SCALE, + $ SUM ) + 310 CONTINUE + ELSE + SCALE = ZERO + SUM = ONE + DO 320 J = 1, N + CALL DLASSQ( M-J+1, A( J, J ), 1, SCALE, SUM ) + 320 CONTINUE + END IF + END IF + VALUE = SCALE*SQRT( SUM ) + END IF +* + DLANTR = VALUE + RETURN +* +* End of DLANTR +* + END