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

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