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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlasr_(char *side, char *pivot, char *direct, integer *m,
   4          integer *n, doublereal *c__, doublereal *s, doublereal *a, integer *
   5         lda)
   6 {
   7     /* System generated locals */
   8     integer a_dim1, a_offset, i__1, i__2;
   9
  10     /* Local variables */
  11     integer i__, j, info;
  12     doublereal temp;
  13     extern logical lsame_(char *, char *);
  14     doublereal ctemp, stemp;
  15     extern /* Subroutine */ int xerbla_(char *, integer *);
  16
  17
  18 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  19 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  20 /*     November 2006 */
  21
  22 /*     .. Scalar Arguments .. */
  23 /*     .. */
  24 /*     .. Array Arguments .. */
  25 /*     .. */
  26
  27 /*  Purpose */
  28 /*  ======= */
  29
  30 /*  DLASR applies a sequence of plane rotations to a real matrix A, */
  31 /*  from either the left or the right. */
  32
  33 /*  When SIDE = 'L', the transformation takes the form */
  34
  35 /*     A := P*A */
  36
  37 /*  and when SIDE = 'R', the transformation takes the form */
  38
  39 /*     A := A*P**T */
  40
  41 /*  where P is an orthogonal matrix consisting of a sequence of z plane */
  42 /*  rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R', */
  43 /*  and P**T is the transpose of P. */
  44
  45 /*  When DIRECT = 'F' (Forward sequence), then */
  46
  47 /*     P = P(z-1) * ... * P(2) * P(1) */
  48
  49 /*  and when DIRECT = 'B' (Backward sequence), then */
  50
  51 /*     P = P(1) * P(2) * ... * P(z-1) */
  52
  53 /*  where P(k) is a plane rotation matrix defined by the 2-by-2 rotation */
  54
  55 /*     R(k) = (  c(k)  s(k) ) */
  56 /*          = ( -s(k)  c(k) ). */
  57
  58 /*  When PIVOT = 'V' (Variable pivot), the rotation is performed */
  59 /*  for the plane (k,k+1), i.e., P(k) has the form */
  60
  61 /*     P(k) = (  1                                            ) */
  62 /*            (       ...                                     ) */
  63 /*            (              1                                ) */
  64 /*            (                   c(k)  s(k)                  ) */
  65 /*            (                  -s(k)  c(k)                  ) */
  66 /*            (                                1              ) */
  67 /*            (                                     ...       ) */
  68 /*            (                                            1  ) */
  69
  70 /*  where R(k) appears as a rank-2 modification to the identity matrix in */
  71 /*  rows and columns k and k+1. */
  72
  73 /*  When PIVOT = 'T' (Top pivot), the rotation is performed for the */
  74 /*  plane (1,k+1), so P(k) has the form */
  75
  76 /*     P(k) = (  c(k)                    s(k)                 ) */
  77 /*            (         1                                     ) */
  78 /*            (              ...                              ) */
  79 /*            (                     1                         ) */
  80 /*            ( -s(k)                    c(k)                 ) */
  81 /*            (                                 1             ) */
  82 /*            (                                      ...      ) */
  83 /*            (                                             1 ) */
  84
  85 /*  where R(k) appears in rows and columns 1 and k+1. */
  86
  87 /*  Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is */
  88 /*  performed for the plane (k,z), giving P(k) the form */
  89
  90 /*     P(k) = ( 1                                             ) */
  91 /*            (      ...                                      ) */
  92 /*            (             1                                 ) */
  93 /*            (                  c(k)                    s(k) ) */
  94 /*            (                         1                     ) */
  95 /*            (                              ...              ) */
  96 /*            (                                     1         ) */
  97 /*            (                 -s(k)                    c(k) ) */
  98
  99 /*  where R(k) appears in rows and columns k and z.  The rotations are */
 100 /*  performed without ever forming P(k) explicitly. */
 101
 102 /*  Arguments */
 103 /*  ========= */
 104
 105 /*  SIDE    (input) CHARACTER*1 */
 106 /*          Specifies whether the plane rotation matrix P is applied to */
 107 /*          A on the left or the right. */
 108 /*          = 'L':  Left, compute A := P*A */
 109 /*          = 'R':  Right, compute A:= A*P**T */
 110
 111 /*  PIVOT   (input) CHARACTER*1 */
 112 /*          Specifies the plane for which P(k) is a plane rotation */
 113 /*          matrix. */
 114 /*          = 'V':  Variable pivot, the plane (k,k+1) */
 115 /*          = 'T':  Top pivot, the plane (1,k+1) */
 116 /*          = 'B':  Bottom pivot, the plane (k,z) */
 117
 118 /*  DIRECT  (input) CHARACTER*1 */
 119 /*          Specifies whether P is a forward or backward sequence of */
 120 /*          plane rotations. */
 121 /*          = 'F':  Forward, P = P(z-1)*...*P(2)*P(1) */
 122 /*          = 'B':  Backward, P = P(1)*P(2)*...*P(z-1) */
 123
 124 /*  M       (input) INTEGER */
 125 /*          The number of rows of the matrix A.  If m <= 1, an immediate */
 126 /*          return is effected. */
 127
 128 /*  N       (input) INTEGER */
 129 /*          The number of columns of the matrix A.  If n <= 1, an */
 130 /*          immediate return is effected. */
 131
 132 /*  C       (input) DOUBLE PRECISION array, dimension */
 133 /*                  (M-1) if SIDE = 'L' */
 134 /*                  (N-1) if SIDE = 'R' */
 135 /*          The cosines c(k) of the plane rotations. */
 136
 137 /*  S       (input) DOUBLE PRECISION array, dimension */
 138 /*                  (M-1) if SIDE = 'L' */
 139 /*                  (N-1) if SIDE = 'R' */
 140 /*          The sines s(k) of the plane rotations.  The 2-by-2 plane */
 141 /*          rotation part of the matrix P(k), R(k), has the form */
 142 /*          R(k) = (  c(k)  s(k) ) */
 143 /*                 ( -s(k)  c(k) ). */
 144
 145 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
 146 /*          The M-by-N matrix A.  On exit, A is overwritten by P*A if */
 147 /*          SIDE = 'R' or by A*P**T if SIDE = 'L'. */
 148
 149 /*  LDA     (input) INTEGER */
 150 /*          The leading dimension of the array A.  LDA >= max(1,M). */
 151
 152 /*  ===================================================================== */
 153
 154 /*     .. Parameters .. */
 155 /*     .. */
 156 /*     .. Local Scalars .. */
 157 /*     .. */
 158 /*     .. External Functions .. */
 159 /*     .. */
 160 /*     .. External Subroutines .. */
 161 /*     .. */
 162 /*     .. Intrinsic Functions .. */
 163 /*     .. */
 164 /*     .. Executable Statements .. */
 165
 166 /*     Test the input parameters */
 167
 168     /* Parameter adjustments */
 169     --c__;
 170     --s;
 171     a_dim1 = *lda;
 172     a_offset = 1 + a_dim1;
 173     a -= a_offset;
 174
 175     /* Function Body */
 176     info = 0;
 177     if (! (lsame_(side, "L") || lsame_(side, "R"))) {
 178         info = 1;
 179     } else if (! (lsame_(pivot, "V") || lsame_(pivot,
 180             "T") || lsame_(pivot, "B"))) {
 181         info = 2;
 182     } else if (! (lsame_(direct, "F") || lsame_(direct,
 183             "B"))) {
 184         info = 3;
 185     } else if (*m < 0) {
 186         info = 4;
 187     } else if (*n < 0) {
 188         info = 5;
 189     } else if (*lda < max(1,*m)) {
 190         info = 9;
 191     }
 192     if (info != 0) {
 193         xerbla_("DLASR ", &info);
 194         return 0;
 195     }
 196
 197 /*     Quick return if possible */
 198
 199     if (*m == 0 || *n == 0) {
 200         return 0;
 201     }
 202     if (lsame_(side, "L")) {
 203
 204 /*        Form  P * A */
 205
 206         if (lsame_(pivot, "V")) {
 207             if (lsame_(direct, "F")) {
 208                 i__1 = *m - 1;
 209                 for (j = 1; j <= i__1; ++j) {
 210                     ctemp = c__[j];
 211                     stemp = s[j];
 212                     if (ctemp != 1. || stemp != 0.) {
 213                         i__2 = *n;
 214                         for (i__ = 1; i__ <= i__2; ++i__) {
 215                             temp = a[j + 1 + i__ * a_dim1];
 216                             a[j + 1 + i__ * a_dim1] = ctemp * temp - stemp *
 217                                     a[j + i__ * a_dim1];
 218                             a[j + i__ * a_dim1] = stemp * temp + ctemp * a[j
 219                                     + i__ * a_dim1];
 220 /* L10: */
 221                         }
 222                     }
 223 /* L20: */
 224                 }
 225             } else if (lsame_(direct, "B")) {
 226                 for (j = *m - 1; j >= 1; --j) {
 227                     ctemp = c__[j];
 228                     stemp = s[j];
 229                     if (ctemp != 1. || stemp != 0.) {
 230                         i__1 = *n;
 231                         for (i__ = 1; i__ <= i__1; ++i__) {
 232                             temp = a[j + 1 + i__ * a_dim1];
 233                             a[j + 1 + i__ * a_dim1] = ctemp * temp - stemp *
 234                                     a[j + i__ * a_dim1];
 235                             a[j + i__ * a_dim1] = stemp * temp + ctemp * a[j
 236                                     + i__ * a_dim1];
 237 /* L30: */
 238                         }
 239                     }
 240 /* L40: */
 241                 }
 242             }
 243         } else if (lsame_(pivot, "T")) {
 244             if (lsame_(direct, "F")) {
 245                 i__1 = *m;
 246                 for (j = 2; j <= i__1; ++j) {
 247                     ctemp = c__[j - 1];
 248                     stemp = s[j - 1];
 249                     if (ctemp != 1. || stemp != 0.) {
 250                         i__2 = *n;
 251                         for (i__ = 1; i__ <= i__2; ++i__) {
 252                             temp = a[j + i__ * a_dim1];
 253                             a[j + i__ * a_dim1] = ctemp * temp - stemp * a[
 254                                     i__ * a_dim1 + 1];
 255                             a[i__ * a_dim1 + 1] = stemp * temp + ctemp * a[
 256                                     i__ * a_dim1 + 1];
 257 /* L50: */
 258                         }
 259                     }
 260 /* L60: */
 261                 }
 262             } else if (lsame_(direct, "B")) {
 263                 for (j = *m; j >= 2; --j) {
 264                     ctemp = c__[j - 1];
 265                     stemp = s[j - 1];
 266                     if (ctemp != 1. || stemp != 0.) {
 267                         i__1 = *n;
 268                         for (i__ = 1; i__ <= i__1; ++i__) {
 269                             temp = a[j + i__ * a_dim1];
 270                             a[j + i__ * a_dim1] = ctemp * temp - stemp * a[
 271                                     i__ * a_dim1 + 1];
 272                             a[i__ * a_dim1 + 1] = stemp * temp + ctemp * a[
 273                                     i__ * a_dim1 + 1];
 274 /* L70: */
 275                         }
 276                     }
 277 /* L80: */
 278                 }
 279             }
 280         } else if (lsame_(pivot, "B")) {
 281             if (lsame_(direct, "F")) {
 282                 i__1 = *m - 1;
 283                 for (j = 1; j <= i__1; ++j) {
 284                     ctemp = c__[j];
 285                     stemp = s[j];
 286                     if (ctemp != 1. || stemp != 0.) {
 287                         i__2 = *n;
 288                         for (i__ = 1; i__ <= i__2; ++i__) {
 289                             temp = a[j + i__ * a_dim1];
 290                             a[j + i__ * a_dim1] = stemp * a[*m + i__ * a_dim1]
 291                                      + ctemp * temp;
 292                             a[*m + i__ * a_dim1] = ctemp * a[*m + i__ *
 293                                     a_dim1] - stemp * temp;
 294 /* L90: */
 295                         }
 296                     }
 297 /* L100: */
 298                 }
 299             } else if (lsame_(direct, "B")) {
 300                 for (j = *m - 1; j >= 1; --j) {
 301                     ctemp = c__[j];
 302                     stemp = s[j];
 303                     if (ctemp != 1. || stemp != 0.) {
 304                         i__1 = *n;
 305                         for (i__ = 1; i__ <= i__1; ++i__) {
 306                             temp = a[j + i__ * a_dim1];
 307                             a[j + i__ * a_dim1] = stemp * a[*m + i__ * a_dim1]
 308                                      + ctemp * temp;
 309                             a[*m + i__ * a_dim1] = ctemp * a[*m + i__ *
 310                                     a_dim1] - stemp * temp;
 311 /* L110: */
 312                         }
 313                     }
 314 /* L120: */
 315                 }
 316             }
 317         }
 318     } else if (lsame_(side, "R")) {
 319
 320 /*        Form A * P' */
 321
 322         if (lsame_(pivot, "V")) {
 323             if (lsame_(direct, "F")) {
 324                 i__1 = *n - 1;
 325                 for (j = 1; j <= i__1; ++j) {
 326                     ctemp = c__[j];
 327                     stemp = s[j];
 328                     if (ctemp != 1. || stemp != 0.) {
 329                         i__2 = *m;
 330                         for (i__ = 1; i__ <= i__2; ++i__) {
 331                             temp = a[i__ + (j + 1) * a_dim1];
 332                             a[i__ + (j + 1) * a_dim1] = ctemp * temp - stemp *
 333                                      a[i__ + j * a_dim1];
 334                             a[i__ + j * a_dim1] = stemp * temp + ctemp * a[
 335                                     i__ + j * a_dim1];
 336 /* L130: */
 337                         }
 338                     }
 339 /* L140: */
 340                 }
 341             } else if (lsame_(direct, "B")) {
 342                 for (j = *n - 1; j >= 1; --j) {
 343                     ctemp = c__[j];
 344                     stemp = s[j];
 345                     if (ctemp != 1. || stemp != 0.) {
 346                         i__1 = *m;
 347                         for (i__ = 1; i__ <= i__1; ++i__) {
 348                             temp = a[i__ + (j + 1) * a_dim1];
 349                             a[i__ + (j + 1) * a_dim1] = ctemp * temp - stemp *
 350                                      a[i__ + j * a_dim1];
 351                             a[i__ + j * a_dim1] = stemp * temp + ctemp * a[
 352                                     i__ + j * a_dim1];
 353 /* L150: */
 354                         }
 355                     }
 356 /* L160: */
 357                 }
 358             }
 359         } else if (lsame_(pivot, "T")) {
 360             if (lsame_(direct, "F")) {
 361                 i__1 = *n;
 362                 for (j = 2; j <= i__1; ++j) {
 363                     ctemp = c__[j - 1];
 364                     stemp = s[j - 1];
 365                     if (ctemp != 1. || stemp != 0.) {
 366                         i__2 = *m;
 367                         for (i__ = 1; i__ <= i__2; ++i__) {
 368                             temp = a[i__ + j * a_dim1];
 369                             a[i__ + j * a_dim1] = ctemp * temp - stemp * a[
 370                                     i__ + a_dim1];
 371                             a[i__ + a_dim1] = stemp * temp + ctemp * a[i__ +
 372                                     a_dim1];
 373 /* L170: */
 374                         }
 375                     }
 376 /* L180: */
 377                 }
 378             } else if (lsame_(direct, "B")) {
 379                 for (j = *n; j >= 2; --j) {
 380                     ctemp = c__[j - 1];
 381                     stemp = s[j - 1];
 382                     if (ctemp != 1. || stemp != 0.) {
 383                         i__1 = *m;
 384                         for (i__ = 1; i__ <= i__1; ++i__) {
 385                             temp = a[i__ + j * a_dim1];
 386                             a[i__ + j * a_dim1] = ctemp * temp - stemp * a[
 387                                     i__ + a_dim1];
 388                             a[i__ + a_dim1] = stemp * temp + ctemp * a[i__ +
 389                                     a_dim1];
 390 /* L190: */
 391                         }
 392                     }
 393 /* L200: */
 394                 }
 395             }
 396         } else if (lsame_(pivot, "B")) {
 397             if (lsame_(direct, "F")) {
 398                 i__1 = *n - 1;
 399                 for (j = 1; j <= i__1; ++j) {
 400                     ctemp = c__[j];
 401                     stemp = s[j];
 402                     if (ctemp != 1. || stemp != 0.) {
 403                         i__2 = *m;
 404                         for (i__ = 1; i__ <= i__2; ++i__) {
 405                             temp = a[i__ + j * a_dim1];
 406                             a[i__ + j * a_dim1] = stemp * a[i__ + *n * a_dim1]
 407                                      + ctemp * temp;
 408                             a[i__ + *n * a_dim1] = ctemp * a[i__ + *n *
 409                                     a_dim1] - stemp * temp;
 410 /* L210: */
 411                         }
 412                     }
 413 /* L220: */
 414                 }
 415             } else if (lsame_(direct, "B")) {
 416                 for (j = *n - 1; j >= 1; --j) {
 417                     ctemp = c__[j];
 418                     stemp = s[j];
 419                     if (ctemp != 1. || stemp != 0.) {
 420                         i__1 = *m;
 421                         for (i__ = 1; i__ <= i__1; ++i__) {
 422                             temp = a[i__ + j * a_dim1];
 423                             a[i__ + j * a_dim1] = stemp * a[i__ + *n * a_dim1]
 424                                      + ctemp * temp;
 425                             a[i__ + *n * a_dim1] = ctemp * a[i__ + *n *
 426                                     a_dim1] - stemp * temp;
 427 /* L230: */
 428                         }
 429                     }
 430 /* L240: */
 431                 }
 432             }
 433         }
 434     }
 435
 436     return 0;
 437
 438 /*     End of DLASR */
 439
 440 } /* dlasr_ */