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

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