src/gmxlib/gmx_lapack/slasdq.c

   1 #include <ctype.h>
   2
   3 #include "gmx_blas.h"
   4 #include "gmx_lapack.h"
   5
   6
   7 void
   8 F77_FUNC(slasdq,SLASDQ)(const char *uplo,
   9                         int *sqre,
  10                         int *n,
  11                         int *ncvt,
  12                         int *nru,
  13                         int *ncc,
  14                         float *d__,
  15                         float *e,
  16                         float *vt,
  17                         int *ldvt,
  18                         float *u,
  19                         int *ldu,
  20                         float *c__,
  21                         int *ldc,
  22                         float *work,
  23                         int *info)
  24 {
  25     const char xuplo=toupper(*uplo);
  26     int c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
  27             i__2;
  28     int c__1 = 1;
  29     int itmp1,itmp2;
  30     int i__, j;
  31     float r__, cs, sn;
  32     int np1, isub;
  33     float smin;
  34     int sqre1;
  35     int iuplo;
  36     int rotate;
  37
  38     --d__;
  39     --e;
  40     vt_dim1 = *ldvt;
  41     vt_offset = 1 + vt_dim1;
  42     vt -= vt_offset;
  43     u_dim1 = *ldu;
  44     u_offset = 1 + u_dim1;
  45     u -= u_offset;
  46     c_dim1 = *ldc;
  47     c_offset = 1 + c_dim1;
  48     c__ -= c_offset;
  49     --work;
  50
  51     *info = 0;
  52     iuplo = 0;
  53     if (xuplo == 'U') {
  54         iuplo = 1;
  55     }
  56     if (xuplo == 'L') {
  57         iuplo = 2;
  58     }
  59
  60     itmp1 = (*n > 1) ? *n : 1;
  61     itmp2 = (*nru > 1) ? *nru : 1;
  62     if (iuplo == 0) {
  63         *info = -1;
  64     } else if (*sqre < 0 || *sqre > 1) {
  65         *info = -2;
  66     } else if (*n < 0) {
  67         *info = -3;
  68     } else if (*ncvt < 0) {
  69         *info = -4;
  70     } else if (*nru < 0) {
  71         *info = -5;
  72     } else if (*ncc < 0) {
  73         *info = -6;
  74     } else if ( (*ncvt == 0 && *ldvt < 1) || (*ncvt > 0 && *ldvt < itmp1)) {
  75         *info = -10;
  76     } else if (*ldu < itmp2) {
  77         *info = -12;
  78     } else if ((*ncc == 0 && *ldc < 1) || (*ncc > 0 && *ldc < itmp1)) {
  79         *info = -14;
  80     }
  81     if (*info != 0) {
  82         return;
  83     }
  84     if (*n == 0) {
  85         return;
  86     }
  87
  88     rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
  89     np1 = *n + 1;
  90     sqre1 = *sqre;
  91
  92     if (iuplo == 1 && sqre1 == 1) {
  93         i__1 = *n - 1;
  94         for (i__ = 1; i__ <= i__1; ++i__) {
  95             F77_FUNC(slartg,SLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
  96             d__[i__] = r__;
  97             e[i__] = sn * d__[i__ + 1];
  98             d__[i__ + 1] = cs * d__[i__ + 1];
  99             if (rotate) {
 100                 work[i__] = cs;
 101                 work[*n + i__] = sn;
 102             }
 103         }
 104         F77_FUNC(slartg,SLARTG)(&d__[*n], &e[*n], &cs, &sn, &r__);
 105         d__[*n] = r__;
 106         e[*n] = 0.f;
 107         if (rotate) {
 108             work[*n] = cs;
 109             work[*n + *n] = sn;
 110         }
 111         iuplo = 2;
 112         sqre1 = 0;
 113
 114         if (*ncvt > 0) {
 115             F77_FUNC(slasr,SLASR)("L", "V", "F", &np1, ncvt, &work[1], &work[np1], &vt[
 116                     vt_offset], ldvt);
 117         }
 118     }
 119     if (iuplo == 2) {
 120         i__1 = *n - 1;
 121         for (i__ = 1; i__ <= i__1; ++i__) {
 122             F77_FUNC(slartg,SLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
 123             d__[i__] = r__;
 124             e[i__] = sn * d__[i__ + 1];
 125             d__[i__ + 1] = cs * d__[i__ + 1];
 126             if (rotate) {
 127                 work[i__] = cs;
 128                 work[*n + i__] = sn;
 129             }
 130         }
 131
 132         if (sqre1 == 1) {
 133             F77_FUNC(slartg,SLARTG)(&d__[*n], &e[*n], &cs, &sn, &r__);
 134             d__[*n] = r__;
 135             if (rotate) {
 136                 work[*n] = cs;
 137                 work[*n + *n] = sn;
 138             }
 139         }
 140         if (*nru > 0) {
 141             if (sqre1 == 0) {
 142                 F77_FUNC(slasr,SLASR)("R", "V", "F", nru, n, &work[1], &work[np1], &u[
 143                         u_offset], ldu);
 144             } else {
 145                 F77_FUNC(slasr,SLASR)("R", "V", "F", nru, &np1, &work[1], &work[np1], &u[
 146                         u_offset], ldu);
 147             }
 148         }
 149         if (*ncc > 0) {
 150             if (sqre1 == 0) {
 151                 F77_FUNC(slasr,SLASR)("L", "V", "F", n, ncc, &work[1], &work[np1], &c__[
 152                         c_offset], ldc);
 153             } else {
 154                 F77_FUNC(slasr,SLASR)("L", "V", "F", &np1, ncc, &work[1], &work[np1], &c__[
 155                         c_offset], ldc);
 156             }
 157         }
 158     }
 159
 160     F77_FUNC(sbdsqr,SBDSQR)("U", n, ncvt, nru, ncc, &d__[1], &e[1], &vt[vt_offset], ldvt, &u[
 161             u_offset], ldu, &c__[c_offset], ldc, &work[1], info);
 162
 163     i__1 = *n;
 164     for (i__ = 1; i__ <= i__1; ++i__) {
 165
 166         isub = i__;
 167         smin = d__[i__];
 168         i__2 = *n;
 169         for (j = i__ + 1; j <= i__2; ++j) {
 170             if (d__[j] < smin) {
 171                 isub = j;
 172                 smin = d__[j];
 173             }
 174         }
 175         if (isub != i__) {
 176             d__[isub] = d__[i__];
 177             d__[i__] = smin;
 178             if (*ncvt > 0) {
 179                 F77_FUNC(sswap,SSWAP)(ncvt, &vt[isub + vt_dim1], ldvt, &vt[i__ + vt_dim1],
 180                         ldvt);
 181             }
 182             if (*nru > 0) {
 183                 F77_FUNC(sswap,SSWAP)(nru, &u[isub * u_dim1 + 1], &c__1, &u[i__ * u_dim1 + 1]
 184                         , &c__1);
 185             }
 186             if (*ncc > 0) {
 187                 F77_FUNC(sswap,SSWAP)(ncc, &c__[isub + c_dim1], ldc, &c__[i__ + c_dim1], ldc)
 188                         ;
 189             }
 190         }
 191     }
 192
 193     return;
 194 }
 195
 196