git.maemo.org Git - opencv/blob - 3rdparty/lapack/dsyr.c

   1
   2 /*  -- translated by f2c (version 19940927).
   3    You must link the resulting object file with the libraries:
   4         -lf2c -lm   (in that order)
   5 */
   6
   7 #include "clapack.h"
   8
   9 /* Subroutine */ int dsyr_(char *uplo, integer *n, doublereal *alpha,
  10         doublereal *x, integer *incx, doublereal *a, integer *lda)
  11 {
  12
  13
  14     /* System generated locals */
  15     integer a_dim1, a_offset, i__1, i__2;
  16
  17     /* Local variables */
  18     static integer info;
  19     static doublereal temp;
  20     static integer i, j;
  21     extern logical lsame_(char *, char *);
  22     static integer ix, jx, kx;
  23     extern /* Subroutine */ int xerbla_(char *, integer *);
  24
  25
  26 /*  Purpose
  27     =======
  28
  29     DSYR   performs the symmetric rank 1 operation
  30
  31        A := alpha*x*x' + A,
  32
  33     where alpha is a real scalar, x is an n element vector and A is an
  34     n by n symmetric matrix.
  35
  36     Parameters
  37     ==========
  38
  39     UPLO   - CHARACTER*1.
  40              On entry, UPLO specifies whether the upper or lower
  41              triangular part of the array A is to be referenced as
  42              follows:
  43
  44                 UPLO = 'U' or 'u'   Only the upper triangular part of A
  45                                     is to be referenced.
  46
  47                 UPLO = 'L' or 'l'   Only the lower triangular part of A
  48                                     is to be referenced.
  49
  50              Unchanged on exit.
  51
  52     N      - INTEGER.
  53              On entry, N specifies the order of the matrix A.
  54              N must be at least zero.
  55              Unchanged on exit.
  56
  57     ALPHA  - DOUBLE PRECISION.
  58              On entry, ALPHA specifies the scalar alpha.
  59              Unchanged on exit.
  60
  61     X      - DOUBLE PRECISION array of dimension at least
  62              ( 1 + ( n - 1 )*abs( INCX ) ).
  63              Before entry, the incremented array X must contain the n
  64              element vector x.
  65              Unchanged on exit.
  66
  67     INCX   - INTEGER.
  68              On entry, INCX specifies the increment for the elements of
  69              X. INCX must not be zero.
  70              Unchanged on exit.
  71
  72     A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  73              Before entry with  UPLO = 'U' or 'u', the leading n by n
  74              upper triangular part of the array A must contain the upper
  75
  76              triangular part of the symmetric matrix and the strictly
  77              lower triangular part of A is not referenced. On exit, the
  78              upper triangular part of the array A is overwritten by the
  79              upper triangular part of the updated matrix.
  80              Before entry with UPLO = 'L' or 'l', the leading n by n
  81              lower triangular part of the array A must contain the lower
  82
  83              triangular part of the symmetric matrix and the strictly
  84              upper triangular part of A is not referenced. On exit, the
  85              lower triangular part of the array A is overwritten by the
  86              lower triangular part of the updated matrix.
  87
  88     LDA    - INTEGER.
  89              On entry, LDA specifies the first dimension of A as declared
  90
  91              in the calling (sub) program. LDA must be at least
  92              max( 1, n ).
  93              Unchanged on exit.
  94
  95
  96     Level 2 Blas routine.
  97
  98     -- Written on 22-October-1986.
  99        Jack Dongarra, Argonne National Lab.
 100        Jeremy Du Croz, Nag Central Office.
 101        Sven Hammarling, Nag Central Office.
 102        Richard Hanson, Sandia National Labs.
 103
 104
 105
 106        Test the input parameters.
 107
 108
 109    Parameter adjustments
 110        Function Body */
 111 #define X(I) x[(I)-1]
 112
 113 #define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]
 114
 115     info = 0;
 116     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
 117         info = 1;
 118     } else if (*n < 0) {
 119         info = 2;
 120     } else if (*incx == 0) {
 121         info = 5;
 122     } else if (*lda < max(1,*n)) {
 123         info = 7;
 124     }
 125     if (info != 0) {
 126         xerbla_("DSYR  ", &info);
 127         return 0;
 128     }
 129
 130 /*     Quick return if possible. */
 131
 132     if (*n == 0 || *alpha == 0.) {
 133         return 0;
 134     }
 135
 136 /*     Set the start point in X if the increment is not unity. */
 137
 138     if (*incx <= 0) {
 139         kx = 1 - (*n - 1) * *incx;
 140     } else if (*incx != 1) {
 141         kx = 1;
 142     }
 143
 144 /*     Start the operations. In this version the elements of A are
 145        accessed sequentially with one pass through the triangular part
 146        of A. */
 147
 148     if (lsame_(uplo, "U")) {
 149
 150 /*        Form  A  when A is stored in upper triangle. */
 151
 152         if (*incx == 1) {
 153             i__1 = *n;
 154             for (j = 1; j <= *n; ++j) {
 155                 if (X(j) != 0.) {
 156                     temp = *alpha * X(j);
 157                     i__2 = j;
 158                     for (i = 1; i <= j; ++i) {
 159                         A(i,j) += X(i) * temp;
 160 /* L10: */
 161                     }
 162                 }
 163 /* L20: */
 164             }
 165         } else {
 166             jx = kx;
 167             i__1 = *n;
 168             for (j = 1; j <= *n; ++j) {
 169                 if (X(jx) != 0.) {
 170                     temp = *alpha * X(jx);
 171                     ix = kx;
 172                     i__2 = j;
 173                     for (i = 1; i <= j; ++i) {
 174                         A(i,j) += X(ix) * temp;
 175                         ix += *incx;
 176 /* L30: */
 177                     }
 178                 }
 179                 jx += *incx;
 180 /* L40: */
 181             }
 182         }
 183     } else {
 184
 185 /*        Form  A  when A is stored in lower triangle. */
 186
 187         if (*incx == 1) {
 188             i__1 = *n;
 189             for (j = 1; j <= *n; ++j) {
 190                 if (X(j) != 0.) {
 191                     temp = *alpha * X(j);
 192                     i__2 = *n;
 193                     for (i = j; i <= *n; ++i) {
 194                         A(i,j) += X(i) * temp;
 195 /* L50: */
 196                     }
 197                 }
 198 /* L60: */
 199             }
 200         } else {
 201             jx = kx;
 202             i__1 = *n;
 203             for (j = 1; j <= *n; ++j) {
 204                 if (X(jx) != 0.) {
 205                     temp = *alpha * X(jx);
 206                     ix = jx;
 207                     i__2 = *n;
 208                     for (i = j; i <= *n; ++i) {
 209                         A(i,j) += X(ix) * temp;
 210                         ix += *incx;
 211 /* L70: */
 212                     }
 213                 }
 214                 jx += *incx;
 215 /* L80: */
 216             }
 217         }
 218     }
 219
 220     return 0;
 221
 222 /*     End of DSYR  . */
 223
 224 } /* dsyr_ */
 225