src/gromacs/linearalgebra/gmx_lapack/sgesdd.c

   1 #include <math.h>
   2 #include "gromacs/utility/real.h"
   3
   4
   5 #include "../gmx_blas.h"
   6 #include "../gmx_lapack.h"
   7 #include "lapack_limits.h"
   8
   9 void
  10 F77_FUNC(sgesdd,SGESDD)(const char *jobz,
  11                         int *m,
  12                         int *n,
  13                         float *a,
  14                         int *lda,
  15                         float *s,
  16                         float *u,
  17                         int *ldu,
  18                         float *vt,
  19                         int *ldvt,
  20                         float *work,
  21                         int *lwork,
  22                         int *iwork,
  23                         int *info)
  24 {
  25     int a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
  26
  27     int ie, iu;
  28     float dum[1], eps;
  29     int ivt, iscl;
  30     float anrm;
  31     int idum[1], ierr, itau;
  32     int minmn, wrkbl, itaup, itauq, mnthr;
  33     int wntqa;
  34     int nwork;
  35     int wntqn, wntqo, wntqs;
  36     int bdspac;
  37     float bignum;
  38     int minwrk, ldwrku, maxwrk, ldwkvt;
  39     float smlnum,minval, safemin;
  40     int wntqas, lquery;
  41     int c__0 = 0;
  42     int c__1 = 1;
  43     float zero = 0.0;
  44     float one = 1.0;
  45
  46
  47     a_dim1 = *lda;
  48     a_offset = 1 + a_dim1;
  49     a -= a_offset;
  50     --s;
  51     u_dim1 = *ldu;
  52     u_offset = 1 + u_dim1;
  53     u -= u_offset;
  54     vt_dim1 = *ldvt;
  55     vt_offset = 1 + vt_dim1;
  56     vt -= vt_offset;
  57     --work;
  58     --iwork;
  59
  60     *info = 0;
  61     minmn = (*m < *n) ? *m : *n;
  62     mnthr = (int) (minmn * 11. / 6.);
  63     wntqa  = (*jobz=='a' || *jobz=='A');
  64     wntqs  = (*jobz=='s' || *jobz=='S');
  65     wntqas = wntqa || wntqs;
  66     wntqn = (*jobz=='o' || *jobz=='O');
  67     wntqo = (*jobz=='n' || *jobz=='N');
  68
  69     minwrk = 1;
  70     maxwrk = 1;
  71     lquery = *lwork == -1;
  72
  73     if (*info == 0 && *m > 0 && *n > 0) {
  74         if (*m >= *n) {
  75
  76             if (wntqn) {
  77                 bdspac = *n * 7;
  78             } else {
  79                 bdspac = *n * 3 * *n + (*n << 2);
  80             }
  81             if (*m >= mnthr) {
  82                 if (wntqn) {
  83
  84                     wrkbl = *n * 67;
  85                     i__1 = wrkbl, i__2 = bdspac + *n;
  86                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
  87                     minwrk = bdspac + *n;
  88                 } else {
  89
  90                     wrkbl = *n * 67;
  91                     i__1 = wrkbl, i__2 = *n + (*m << 5);
  92                     wrkbl = (i__1 > i__2) ? i__1 : i__2;
  93                     i__1 = wrkbl, i__2 = bdspac + *n * 3;
  94                     wrkbl = (i__1 > i__2) ? i__1 : i__2;
  95                     maxwrk = wrkbl + *n * *n;
  96                     minwrk = bdspac + *n * *n + *n * 3;
  97                 }
  98             } else {
  99
 100                 wrkbl = *n * 3 + (*m + *n*32);
 101                 if (wntqn) {
 102                     i__1 = wrkbl, i__2 = bdspac + *n * 3;
 103                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 104                     minwrk = *n * 3 + ((*m > bdspac) ? *m : bdspac);
 105                 } else {
 106                     i__1 = maxwrk, i__2 = bdspac + *n * 3;
 107                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 108                     minwrk = *n * 3 + ((*m > bdspac) ? *m : bdspac);
 109                 }
 110             }
 111         } else {
 112
 113             if (wntqn) {
 114                 bdspac = *m * 7;
 115             } else {
 116                 bdspac = *m * 3 * *m + (*m*4);
 117             }
 118             if (*n >= mnthr) {
 119                 if (wntqn) {
 120
 121                     wrkbl = *m * 67;
 122                     i__1 = wrkbl, i__2 = bdspac + *m;
 123                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 124                     minwrk = bdspac + *m;
 125                 } else {
 126
 127                     wrkbl = *m * 67;
 128                     i__1 = wrkbl, i__2 = *m + (*n*32);
 129                     wrkbl = (i__1 > i__2) ? i__1 : i__2;
 130
 131                     i__1 = wrkbl, i__2 = bdspac + *m * 3;
 132                     wrkbl = (i__1 > i__2) ? i__1 : i__2;
 133                     maxwrk = wrkbl + *m * *m;
 134                     minwrk = bdspac + *m * *m + *m * 3;
 135                 }
 136             } else {
 137                 wrkbl = *m * 3 + (*m + *n*32);
 138                 if (wntqn) {
 139                     i__1 = wrkbl, i__2 = bdspac + *m * 3;
 140                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 141                     minwrk = *m * 3 + ((*m > bdspac) ? *m : bdspac);
 142                 } else {
 143                     i__1 = wrkbl, i__2 = bdspac + *m * 3;
 144                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 145                     minwrk = *m * 3 + ((*m > bdspac) ? *m : bdspac);
 146                 }
 147             }
 148         }
 149         work[1] = (float) maxwrk;
 150     }
 151
 152     if( lquery != 0)
 153     {
 154         return;
 155     }
 156
 157     if (*m == 0 || *n == 0) {
 158         if (*lwork >= 1) {
 159             work[1] = 1.;
 160         }
 161         return;
 162     }
 163     eps = GMX_FLOAT_EPS;
 164     minval = GMX_FLOAT_MIN;
 165     safemin = minval / eps;
 166     smlnum = sqrt(safemin) / eps;
 167
 168
 169     bignum = 1. / smlnum;
 170
 171
 172     anrm = F77_FUNC(slange,SLANGE)("M", m, n, &a[a_offset], lda, dum);
 173     iscl = 0;
 174     if (anrm > 0. && anrm < smlnum) {
 175         iscl = 1;
 176         F77_FUNC(slascl,SLASCL)("G",&c__0,&c__0,&anrm,&smlnum,m,n,&a[a_offset],lda,&ierr);
 177     } else if (anrm > bignum) {
 178         iscl = 1;
 179         F77_FUNC(slascl,SLASCL)("G",&c__0,&c__0,&anrm,&bignum,m,n,&a[a_offset],lda,&ierr);
 180     }
 181
 182     if (*m >= *n) {
 183         if (*m >= mnthr) {
 184
 185             if (wntqn) {
 186
 187                 itau = 1;
 188                 nwork = itau + *n;
 189
 190                 i__1 = *lwork - nwork + 1;
 191                 F77_FUNC(sgeqrf,SGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 192                         i__1, &ierr);
 193
 194                 i__1 = *n - 1;
 195                 i__2 = *n - 1;
 196                 F77_FUNC(slaset,SLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
 197                         lda);
 198                 ie = 1;
 199                 itauq = ie + *n;
 200                 itaup = itauq + *n;
 201                 nwork = itaup + *n;
 202
 203                 i__1 = *lwork - nwork + 1;
 204                 F77_FUNC(sgebrd,SGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
 205                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 206                 nwork = ie + *n;
 207
 208                 F77_FUNC(sbdsdc,SBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 209                          dum, idum, &work[nwork], &iwork[1], info);
 210
 211             } else {
 212                 iu = 1;
 213
 214                 ldwrku = *n;
 215                 itau = iu + ldwrku * *n;
 216                 nwork = itau + *n;
 217
 218                 i__1 = *lwork - nwork + 1;
 219                 F77_FUNC(sgeqrf,SGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 220                         i__1, &ierr);
 221                 F77_FUNC(slacpy,SLACPY)("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
 222
 223                 i__1 = *lwork - nwork + 1;
 224                 F77_FUNC(sorgqr,SORGQR)(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
 225                          &i__1, &ierr);
 226
 227                 i__1 = *n - 1;
 228                 i__2 = *n - 1;
 229                 F77_FUNC(slaset,SLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
 230                         lda);
 231                 ie = itau;
 232                 itauq = ie + *n;
 233                 itaup = itauq + *n;
 234                 nwork = itaup + *n;
 235
 236                 i__1 = *lwork - nwork + 1;
 237                 F77_FUNC(sgebrd,SGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
 238                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 239
 240                 F77_FUNC(sbdsdc,SBDSDC)("U", "I", n, &s[1], &work[ie], &work[iu], n, &vt[
 241                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 242                         info);
 243
 244                 i__1 = *lwork - nwork + 1;
 245                 F77_FUNC(sormbr,SORMBR)("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
 246                         itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
 247                         ierr);
 248                 i__1 = *lwork - nwork + 1;
 249                 F77_FUNC(sormbr,SORMBR)("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
 250                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
 251                         ierr);
 252
 253                 F77_FUNC(sgemm,SGEMM)("N", "N", m, n, n, &one, &u[u_offset], ldu, &work[iu]
 254                         , &ldwrku, &zero, &a[a_offset], lda);
 255
 256                 F77_FUNC(slacpy,SLACPY)("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
 257
 258             }
 259
 260         } else {
 261             ie = 1;
 262             itauq = ie + *n;
 263             itaup = itauq + *n;
 264             nwork = itaup + *n;
 265
 266             i__1 = *lwork - nwork + 1;
 267             F77_FUNC(sgebrd,SGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
 268                     work[itaup], &work[nwork], &i__1, &ierr);
 269             if (wntqn) {
 270
 271                 F77_FUNC(sbdsdc,SBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 272                          dum, idum, &work[nwork], &iwork[1], info);
 273             } else {
 274
 275                 F77_FUNC(slaset,SLASET)("F", m, m, &zero, &zero, &u[u_offset], ldu);
 276                 F77_FUNC(sbdsdc,SBDSDC)("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
 277                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 278                         info);
 279
 280                 i__1 = *m - *n;
 281                 i__2 = *m - *n;
 282                 F77_FUNC(slaset,SLASET)("F", &i__1, &i__2, &zero, &one, &u[*n + 1 + (*n +
 283                         1) * u_dim1], ldu);
 284
 285                 i__1 = *lwork - nwork + 1;
 286                 F77_FUNC(sormbr,SORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
 287                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 288                 i__1 = *lwork - nwork + 1;
 289                 F77_FUNC(sormbr,SORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
 290                         itaup], &vt[vt_offset],ldvt,&work[nwork],&i__1,&ierr);
 291             }
 292
 293         }
 294
 295     } else {
 296
 297         if (*n >= mnthr) {
 298
 299             if (wntqn) {
 300
 301                 itau = 1;
 302                 nwork = itau + *m;
 303
 304                 i__1 = *lwork - nwork + 1;
 305                 F77_FUNC(sgelqf,SGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 306                         i__1, &ierr);
 307
 308                 i__1 = *m - 1;
 309                 i__2 = *m - 1;
 310                 F77_FUNC(slaset,SLASET)("U", &i__1, &i__2, &zero, &zero, &a[(a_dim1*2) +
 311                         1], lda);
 312                 ie = 1;
 313                 itauq = ie + *m;
 314                 itaup = itauq + *m;
 315                 nwork = itaup + *m;
 316
 317                 i__1 = *lwork - nwork + 1;
 318                 F77_FUNC(sgebrd,SGEBRD)(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
 319                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 320                 nwork = ie + *m;
 321
 322                 F77_FUNC(sbdsdc,SBDSDC)("U", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 323                          dum, idum, &work[nwork], &iwork[1], info);
 324
 325             } else {
 326
 327                 ivt = 1;
 328
 329                 ldwkvt = *m;
 330                 itau = ivt + ldwkvt * *m;
 331                 nwork = itau + *m;
 332
 333                 i__1 = *lwork - nwork + 1;
 334                 F77_FUNC(sgelqf,SGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 335                         i__1, &ierr);
 336                 F77_FUNC(slacpy,SLACPY)("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
 337
 338                 i__1 = *lwork - nwork + 1;
 339                 F77_FUNC(sorglq,SORGLQ)(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
 340                         nwork], &i__1, &ierr);
 341
 342                 i__1 = *m - 1;
 343                 i__2 = *m - 1;
 344                 F77_FUNC(slaset,SLASET)("U", &i__1, &i__2, &zero, &zero, &a[(a_dim1*2) +
 345                         1], lda);
 346                 ie = itau;
 347                 itauq = ie + *m;
 348                 itaup = itauq + *m;
 349                 nwork = itaup + *m;
 350
 351                 i__1 = *lwork - nwork + 1;
 352                 F77_FUNC(sgebrd,SGEBRD)(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
 353                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 354
 355                 F77_FUNC(sbdsdc,SBDSDC)("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
 356                         work[ivt], &ldwkvt, dum, idum, &work[nwork], &iwork[1]
 357                         , info);
 358
 359                 i__1 = *lwork - nwork + 1;
 360                 F77_FUNC(sormbr,SORMBR)("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
 361                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 362                 i__1 = *lwork - nwork + 1;
 363                 F77_FUNC(sormbr,SORMBR)("P", "R", "T", m, m, m, &a[a_offset], lda, &work[
 364                         itaup], &work[ivt], &ldwkvt, &work[nwork], &i__1, &
 365                         ierr);
 366
 367                 F77_FUNC(sgemm,SGEMM)("N", "N", m, n, m, &one, &work[ivt], &ldwkvt, &vt[
 368                         vt_offset], ldvt, &zero, &a[a_offset], lda);
 369
 370                 F77_FUNC(slacpy,SLACPY)("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
 371
 372             }
 373
 374         } else {
 375
 376             ie = 1;
 377             itauq = ie + *m;
 378             itaup = itauq + *m;
 379             nwork = itaup + *m;
 380
 381             i__1 = *lwork - nwork + 1;
 382             F77_FUNC(sgebrd,SGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
 383                     work[itaup], &work[nwork], &i__1, &ierr);
 384             if (wntqn) {
 385
 386                 F77_FUNC(sbdsdc,SBDSDC)("L", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 387                          dum, idum, &work[nwork], &iwork[1], info);
 388             } else {
 389                 F77_FUNC(slaset,SLASET)("F", n, n, &zero, &zero, &vt[vt_offset], ldvt);
 390                 F77_FUNC(sbdsdc,SBDSDC)("L", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &vt[
 391                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 392                         info);
 393
 394                 i__1 = *n - *m;
 395                 i__2 = *n - *m;
 396                 F77_FUNC(slaset,SLASET)("F", &i__1, &i__2, &zero, &one, &vt[*m + 1 + (*m +
 397                         1) * vt_dim1], ldvt);
 398
 399                 i__1 = *lwork - nwork + 1;
 400                 F77_FUNC(sormbr,SORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
 401                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 402                 i__1 = *lwork - nwork + 1;
 403                 F77_FUNC(sormbr,SORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
 404                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
 405                         ierr);
 406             }
 407
 408         }
 409
 410     }
 411
 412     if (iscl == 1) {
 413         if (anrm > bignum) {
 414             F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
 415                     minmn, &ierr);
 416         }
 417         if (anrm < smlnum) {
 418             F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
 419                     minmn, &ierr);
 420         }
 421     }
 422
 423     work[1] = (float) maxwrk;
 424
 425     return;
 426
 427 }
 428
 429