dlavsy.c

Go to the documentation of this file.

/* dlavsy.f -- translated by f2c (version 20061008).
   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"
#include "blaswrap.h"

/* Table of constant values */

static doublereal c_b15 = 1.;
static integer c__1 = 1;

/* Subroutine */ int dlavsy_(char *uplo, char *trans, char *diag, integer *n, 
        integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal 
        *b, integer *ldb, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, i__1;

    /* Local variables */
    integer j, k;
    doublereal t1, t2, d11, d12, d21, d22;
    integer kp;
    extern /* Subroutine */ int dger_(integer *, integer *, doublereal *, 
            doublereal *, integer *, doublereal *, integer *, doublereal *, 
            integer *), dscal_(integer *, doublereal *, doublereal *, integer 
            *);
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int dgemv_(char *, integer *, integer *, 
            doublereal *, doublereal *, integer *, doublereal *, integer *, 
            doublereal *, doublereal *, integer *), dswap_(integer *, 
            doublereal *, integer *, doublereal *, integer *), xerbla_(char *, 
             integer *);
    logical nounit;


/*  -- LAPACK auxiliary routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  DLAVSY  performs one of the matrix-vector operations */
/*     x := A*x  or  x := A'*x, */
/*  where x is an N element vector and A is one of the factors */
/*  from the block U*D*U' or L*D*L' factorization computed by DSYTRF. */

/*  If TRANS = 'N', multiplies by U  or U * D  (or L  or L * D) */
/*  If TRANS = 'T', multiplies by U' or D * U' (or L' or D * L') */
/*  If TRANS = 'C', multiplies by U' or D * U' (or L' or D * L') */

/*  Arguments */
/*  ========= */

/*  UPLO    (input) CHARACTER*1 */
/*          Specifies whether the factor stored in A is upper or lower */
/*          triangular. */
/*          = 'U':  Upper triangular */
/*          = 'L':  Lower triangular */

/*  TRANS   (input) CHARACTER*1 */
/*          Specifies the operation to be performed: */
/*          = 'N':  x := A*x */
/*          = 'T':  x := A'*x */
/*          = 'C':  x := A'*x */

/*  DIAG    (input) CHARACTER*1 */
/*          Specifies whether or not the diagonal blocks are unit */
/*          matrices.  If the diagonal blocks are assumed to be unit, */
/*          then A = U or A = L, otherwise A = U*D or A = L*D. */
/*          = 'U':  Diagonal blocks are assumed to be unit matrices. */
/*          = 'N':  Diagonal blocks are assumed to be non-unit matrices. */

/*  N       (input) INTEGER */
/*          The number of rows and columns of the matrix A.  N >= 0. */

/*  NRHS    (input) INTEGER */
/*          The number of right hand sides, i.e., the number of vectors */
/*          x to be multiplied by A.  NRHS >= 0. */

/*  A       (input) DOUBLE PRECISION array, dimension (LDA,N) */
/*          The block diagonal matrix D and the multipliers used to */
/*          obtain the factor U or L as computed by DSYTRF. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A.  LDA >= max(1,N). */

/*  IPIV    (input) INTEGER array, dimension (N) */
/*          The pivot indices from DSYTRF. */

/*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
/*          On entry, B contains NRHS vectors of length N. */
/*          On exit, B is overwritten with the product A * B. */

/*  LDB     (input) INTEGER */
/*          The leading dimension of the array B.  LDB >= max(1,N). */

/*  INFO    (output) INTEGER */
/*          = 0: successful exit */
/*          < 0: if INFO = -k, the k-th argument had an illegal value */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Test the input parameters. */

    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --ipiv;
    b_dim1 = *ldb;
    b_offset = 1 + b_dim1;
    b -= b_offset;

    /* Function Body */
    *info = 0;
    if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
        *info = -1;
    } else if (! lsame_(trans, "N") && ! lsame_(trans, 
            "T") && ! lsame_(trans, "C")) {
        *info = -2;
    } else if (! lsame_(diag, "U") && ! lsame_(diag, 
            "N")) {
        *info = -3;
    } else if (*n < 0) {
        *info = -4;
    } else if (*lda < max(1,*n)) {
        *info = -6;
    } else if (*ldb < max(1,*n)) {
        *info = -9;
    }
    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("DLAVSY ", &i__1);
        return 0;
    }

/*     Quick return if possible. */

    if (*n == 0) {
        return 0;
    }

    nounit = lsame_(diag, "N");
/* ------------------------------------------ */

/*     Compute  B := A * B  (No transpose) */

/* ------------------------------------------ */
    if (lsame_(trans, "N")) {

/*        Compute  B := U*B */
/*        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) */

        if (lsame_(uplo, "U")) {

/*        Loop forward applying the transformations. */

            k = 1;
L10:
            if (k > *n) {
                goto L30;
            }
            if (ipiv[k] > 0) {

/*              1 x 1 pivot block */

/*              Multiply by the diagonal element if forming U * D. */

                if (nounit) {
                    dscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
                }

/*              Multiply by  P(K) * inv(U(K))  if K > 1. */

                if (k > 1) {

/*                 Apply the transformation. */

                    i__1 = k - 1;
                    dger_(&i__1, nrhs, &c_b15, &a[k * a_dim1 + 1], &c__1, &b[
                            k + b_dim1], ldb, &b[b_dim1 + 1], ldb);

/*                 Interchange if P(K) .ne. I. */

                    kp = ipiv[k];
                    if (kp != k) {
                        dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
                                ldb);
                    }
                }
                ++k;
            } else {

/*              2 x 2 pivot block */

/*              Multiply by the diagonal block if forming U * D. */

                if (nounit) {
                    d11 = a[k + k * a_dim1];
                    d22 = a[k + 1 + (k + 1) * a_dim1];
                    d12 = a[k + (k + 1) * a_dim1];
                    d21 = d12;
                    i__1 = *nrhs;
                    for (j = 1; j <= i__1; ++j) {
                        t1 = b[k + j * b_dim1];
                        t2 = b[k + 1 + j * b_dim1];
                        b[k + j * b_dim1] = d11 * t1 + d12 * t2;
                        b[k + 1 + j * b_dim1] = d21 * t1 + d22 * t2;
/* L20: */
                    }
                }

/*              Multiply by  P(K) * inv(U(K))  if K > 1. */

                if (k > 1) {

/*                 Apply the transformations. */

                    i__1 = k - 1;
                    dger_(&i__1, nrhs, &c_b15, &a[k * a_dim1 + 1], &c__1, &b[
                            k + b_dim1], ldb, &b[b_dim1 + 1], ldb);
                    i__1 = k - 1;
                    dger_(&i__1, nrhs, &c_b15, &a[(k + 1) * a_dim1 + 1], &
                            c__1, &b[k + 1 + b_dim1], ldb, &b[b_dim1 + 1], 
                            ldb);

/*                 Interchange if P(K) .ne. I. */

                    kp = (i__1 = ipiv[k], abs(i__1));
                    if (kp != k) {
                        dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
                                ldb);
                    }
                }
                k += 2;
            }
            goto L10;
L30:

/*        Compute  B := L*B */
/*        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) . */

            ;
        } else {

/*           Loop backward applying the transformations to B. */

            k = *n;
L40:
            if (k < 1) {
                goto L60;
            }

/*           Test the pivot index.  If greater than zero, a 1 x 1 */
/*           pivot was used, otherwise a 2 x 2 pivot was used. */

            if (ipiv[k] > 0) {

/*              1 x 1 pivot block: */

/*              Multiply by the diagonal element if forming L * D. */

                if (nounit) {
                    dscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
                }

/*              Multiply by  P(K) * inv(L(K))  if K < N. */

                if (k != *n) {
                    kp = ipiv[k];

/*                 Apply the transformation. */

                    i__1 = *n - k;
                    dger_(&i__1, nrhs, &c_b15, &a[k + 1 + k * a_dim1], &c__1, 
                            &b[k + b_dim1], ldb, &b[k + 1 + b_dim1], ldb);

/*                 Interchange if a permutation was applied at the */
/*                 K-th step of the factorization. */

                    if (kp != k) {
                        dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
                                ldb);
                    }
                }
                --k;

            } else {

/*              2 x 2 pivot block: */

/*              Multiply by the diagonal block if forming L * D. */

                if (nounit) {
                    d11 = a[k - 1 + (k - 1) * a_dim1];
                    d22 = a[k + k * a_dim1];
                    d21 = a[k + (k - 1) * a_dim1];
                    d12 = d21;
                    i__1 = *nrhs;
                    for (j = 1; j <= i__1; ++j) {
                        t1 = b[k - 1 + j * b_dim1];
                        t2 = b[k + j * b_dim1];
                        b[k - 1 + j * b_dim1] = d11 * t1 + d12 * t2;
                        b[k + j * b_dim1] = d21 * t1 + d22 * t2;
/* L50: */
                    }
                }

/*              Multiply by  P(K) * inv(L(K))  if K < N. */

                if (k != *n) {

/*                 Apply the transformation. */

                    i__1 = *n - k;
                    dger_(&i__1, nrhs, &c_b15, &a[k + 1 + k * a_dim1], &c__1, 
                            &b[k + b_dim1], ldb, &b[k + 1 + b_dim1], ldb);
                    i__1 = *n - k;
                    dger_(&i__1, nrhs, &c_b15, &a[k + 1 + (k - 1) * a_dim1], &
                            c__1, &b[k - 1 + b_dim1], ldb, &b[k + 1 + b_dim1], 
                             ldb);

/*                 Interchange if a permutation was applied at the */
/*                 K-th step of the factorization. */

                    kp = (i__1 = ipiv[k], abs(i__1));
                    if (kp != k) {
                        dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
                                ldb);
                    }
                }
                k += -2;
            }
            goto L40;
L60:
            ;
        }
/* ---------------------------------------- */

/*     Compute  B := A' * B  (transpose) */

/* ---------------------------------------- */
    } else {

/*        Form  B := U'*B */
/*        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1)) */
/*        and   U' = inv(U'(1))*P(1)* ... *inv(U'(m))*P(m) */

        if (lsame_(uplo, "U")) {

/*           Loop backward applying the transformations. */

            k = *n;
L70:
            if (k < 1) {
                goto L90;
            }

/*           1 x 1 pivot block. */

            if (ipiv[k] > 0) {
                if (k > 1) {

/*                 Interchange if P(K) .ne. I. */

                    kp = ipiv[k];
                    if (kp != k) {
                        dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
                                ldb);
                    }

/*                 Apply the transformation */

                    i__1 = k - 1;
                    dgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset], 
                            ldb, &a[k * a_dim1 + 1], &c__1, &c_b15, &b[k + 
                            b_dim1], ldb);
                }
                if (nounit) {
                    dscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
                }
                --k;

/*           2 x 2 pivot block. */

            } else {
                if (k > 2) {

/*                 Interchange if P(K) .ne. I. */

                    kp = (i__1 = ipiv[k], abs(i__1));
                    if (kp != k - 1) {
                        dswap_(nrhs, &b[k - 1 + b_dim1], ldb, &b[kp + b_dim1], 
                                 ldb);
                    }

/*                 Apply the transformations */

                    i__1 = k - 2;
                    dgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset], 
                            ldb, &a[k * a_dim1 + 1], &c__1, &c_b15, &b[k + 
                            b_dim1], ldb);
                    i__1 = k - 2;
                    dgemv_("Transpose", &i__1, nrhs, &c_b15, &b[b_offset], 
                            ldb, &a[(k - 1) * a_dim1 + 1], &c__1, &c_b15, &b[
                            k - 1 + b_dim1], ldb);
                }

/*              Multiply by the diagonal block if non-unit. */

                if (nounit) {
                    d11 = a[k - 1 + (k - 1) * a_dim1];
                    d22 = a[k + k * a_dim1];
                    d12 = a[k - 1 + k * a_dim1];
                    d21 = d12;
                    i__1 = *nrhs;
                    for (j = 1; j <= i__1; ++j) {
                        t1 = b[k - 1 + j * b_dim1];
                        t2 = b[k + j * b_dim1];
                        b[k - 1 + j * b_dim1] = d11 * t1 + d12 * t2;
                        b[k + j * b_dim1] = d21 * t1 + d22 * t2;
/* L80: */
                    }
                }
                k += -2;
            }
            goto L70;
L90:

/*        Form  B := L'*B */
/*        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) */
/*        and   L' = inv(L'(m))*P(m)* ... *inv(L'(1))*P(1) */

            ;
        } else {

/*           Loop forward applying the L-transformations. */

            k = 1;
L100:
            if (k > *n) {
                goto L120;
            }

/*           1 x 1 pivot block */

            if (ipiv[k] > 0) {
                if (k < *n) {

/*                 Interchange if P(K) .ne. I. */

                    kp = ipiv[k];
                    if (kp != k) {
                        dswap_(nrhs, &b[k + b_dim1], ldb, &b[kp + b_dim1], 
                                ldb);
                    }

/*                 Apply the transformation */

                    i__1 = *n - k;
                    dgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 1 + 
                            b_dim1], ldb, &a[k + 1 + k * a_dim1], &c__1, &
                            c_b15, &b[k + b_dim1], ldb);
                }
                if (nounit) {
                    dscal_(nrhs, &a[k + k * a_dim1], &b[k + b_dim1], ldb);
                }
                ++k;

/*           2 x 2 pivot block. */

            } else {
                if (k < *n - 1) {

/*              Interchange if P(K) .ne. I. */

                    kp = (i__1 = ipiv[k], abs(i__1));
                    if (kp != k + 1) {
                        dswap_(nrhs, &b[k + 1 + b_dim1], ldb, &b[kp + b_dim1], 
                                 ldb);
                    }

/*                 Apply the transformation */

                    i__1 = *n - k - 1;
                    dgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 2 + 
                            b_dim1], ldb, &a[k + 2 + (k + 1) * a_dim1], &c__1, 
                             &c_b15, &b[k + 1 + b_dim1], ldb);
                    i__1 = *n - k - 1;
                    dgemv_("Transpose", &i__1, nrhs, &c_b15, &b[k + 2 + 
                            b_dim1], ldb, &a[k + 2 + k * a_dim1], &c__1, &
                            c_b15, &b[k + b_dim1], ldb);
                }

/*              Multiply by the diagonal block if non-unit. */

                if (nounit) {
                    d11 = a[k + k * a_dim1];
                    d22 = a[k + 1 + (k + 1) * a_dim1];
                    d21 = a[k + 1 + k * a_dim1];
                    d12 = d21;
                    i__1 = *nrhs;
                    for (j = 1; j <= i__1; ++j) {
                        t1 = b[k + j * b_dim1];
                        t2 = b[k + 1 + j * b_dim1];
                        b[k + j * b_dim1] = d11 * t1 + d12 * t2;
                        b[k + 1 + j * b_dim1] = d21 * t1 + d22 * t2;
/* L110: */
                    }
                }
                k += 2;
            }
            goto L100;
L120:
            ;
        }

    }
    return 0;

/*     End of DLAVSY */

} /* dlavsy_ */