src/gmxlib/gmx_lapack/sgesdd.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <math.h>
  36 #include <types/simple.h>
  37
  38
  39 #include "gmx_blas.h"
  40 #include "gmx_lapack.h"
  41 #include "lapack_limits.h"
  42
  43 void
  44 F77_FUNC(sgesdd,SGESDD)(const char *jobz,
  45                         int *m,
  46                         int *n,
  47                         float *a,
  48                         int *lda,
  49                         float *s,
  50                         float *u,
  51                         int *ldu,
  52                         float *vt,
  53                         int *ldvt,
  54                         float *work,
  55                         int *lwork,
  56                         int *iwork,
  57                         int *info)
  58 {
  59     int a_dim1, a_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
  60
  61     int ie, iu;
  62     float dum[1], eps;
  63     int ivt, iscl;
  64     float anrm;
  65     int idum[1], ierr, itau;
  66     int minmn, wrkbl, itaup, itauq, mnthr;
  67     int wntqa;
  68     int nwork;
  69     int wntqn, wntqo, wntqs;
  70     int bdspac;
  71     float bignum;
  72     int minwrk, ldwrku, maxwrk, ldwkvt;
  73     float smlnum,minval, safemin;
  74     int wntqas, lquery;
  75     int c__0 = 0;
  76     int c__1 = 1;
  77     float zero = 0.0;
  78     float one = 1.0;
  79
  80
  81     a_dim1 = *lda;
  82     a_offset = 1 + a_dim1;
  83     a -= a_offset;
  84     --s;
  85     u_dim1 = *ldu;
  86     u_offset = 1 + u_dim1;
  87     u -= u_offset;
  88     vt_dim1 = *ldvt;
  89     vt_offset = 1 + vt_dim1;
  90     vt -= vt_offset;
  91     --work;
  92     --iwork;
  93
  94     *info = 0;
  95     minmn = (*m < *n) ? *m : *n;
  96     mnthr = (int) (minmn * 11. / 6.);
  97     wntqa  = (*jobz=='a' || *jobz=='A');
  98     wntqs  = (*jobz=='s' || *jobz=='S');
  99     wntqas = wntqa || wntqs;
 100     wntqn = (*jobz=='o' || *jobz=='O');
 101     wntqo = (*jobz=='n' || *jobz=='N');
 102
 103     minwrk = 1;
 104     maxwrk = 1;
 105     lquery = *lwork == -1;
 106
 107     if (*info == 0 && *m > 0 && *n > 0) {
 108         if (*m >= *n) {
 109
 110             if (wntqn) {
 111                 bdspac = *n * 7;
 112             } else {
 113                 bdspac = *n * 3 * *n + (*n << 2);
 114             }
 115             if (*m >= mnthr) {
 116                 if (wntqn) {
 117
 118                     wrkbl = *n * 67;
 119                     i__1 = wrkbl, i__2 = bdspac + *n;
 120                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 121                     minwrk = bdspac + *n;
 122                 } else {
 123
 124                     wrkbl = *n * 67;
 125                     i__1 = wrkbl, i__2 = *n + (*m << 5);
 126                     wrkbl = (i__1 > i__2) ? i__1 : i__2;
 127                     i__1 = wrkbl, i__2 = bdspac + *n * 3;
 128                     wrkbl = (i__1 > i__2) ? i__1 : i__2;
 129                     maxwrk = wrkbl + *n * *n;
 130                     minwrk = bdspac + *n * *n + *n * 3;
 131                 }
 132             } else {
 133
 134                 wrkbl = *n * 3 + (*m + *n*32);
 135                 if (wntqn) {
 136                     i__1 = wrkbl, i__2 = bdspac + *n * 3;
 137                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 138                     minwrk = *n * 3 + ((*m > bdspac) ? *m : bdspac);
 139                 } else {
 140                     i__1 = maxwrk, i__2 = bdspac + *n * 3;
 141                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 142                     minwrk = *n * 3 + ((*m > bdspac) ? *m : bdspac);
 143                 }
 144             }
 145         } else {
 146
 147             if (wntqn) {
 148                 bdspac = *m * 7;
 149             } else {
 150                 bdspac = *m * 3 * *m + (*m*4);
 151             }
 152             if (*n >= mnthr) {
 153                 if (wntqn) {
 154
 155                     wrkbl = *m * 67;
 156                     i__1 = wrkbl, i__2 = bdspac + *m;
 157                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 158                     minwrk = bdspac + *m;
 159                 } else {
 160
 161                     wrkbl = *m * 67;
 162                     i__1 = wrkbl, i__2 = *m + (*n*32);
 163                     wrkbl = (i__1 > i__2) ? i__1 : i__2;
 164
 165                     i__1 = wrkbl, i__2 = bdspac + *m * 3;
 166                     wrkbl = (i__1 > i__2) ? i__1 : i__2;
 167                     maxwrk = wrkbl + *m * *m;
 168                     minwrk = bdspac + *m * *m + *m * 3;
 169                 }
 170             } else {
 171                 wrkbl = *m * 3 + (*m + *n*32);
 172                 if (wntqn) {
 173                     i__1 = wrkbl, i__2 = bdspac + *m * 3;
 174                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 175                     minwrk = *m * 3 + ((*m > bdspac) ? *m : bdspac);
 176                 } else {
 177                     i__1 = wrkbl, i__2 = bdspac + *m * 3;
 178                     maxwrk = (i__1 > i__2) ? i__1 : i__2;
 179                     minwrk = *m * 3 + ((*m > bdspac) ? *m : bdspac);
 180                 }
 181             }
 182         }
 183         work[1] = (float) maxwrk;
 184     }
 185
 186     if( lquery != 0)
 187     {
 188         return;
 189     }
 190
 191     if (*m == 0 || *n == 0) {
 192         if (*lwork >= 1) {
 193             work[1] = 1.;
 194         }
 195         return;
 196     }
 197     eps = GMX_FLOAT_EPS;
 198     minval = GMX_FLOAT_MIN;
 199     safemin = minval / eps;
 200     smlnum = sqrt(safemin) / eps;
 201
 202
 203     bignum = 1. / smlnum;
 204
 205
 206     anrm = F77_FUNC(slange,SLANGE)("M", m, n, &a[a_offset], lda, dum);
 207     iscl = 0;
 208     if (anrm > 0. && anrm < smlnum) {
 209         iscl = 1;
 210         F77_FUNC(slascl,SLASCL)("G",&c__0,&c__0,&anrm,&smlnum,m,n,&a[a_offset],lda,&ierr);
 211     } else if (anrm > bignum) {
 212         iscl = 1;
 213         F77_FUNC(slascl,SLASCL)("G",&c__0,&c__0,&anrm,&bignum,m,n,&a[a_offset],lda,&ierr);
 214     }
 215
 216     if (*m >= *n) {
 217         if (*m >= mnthr) {
 218
 219             if (wntqn) {
 220
 221                 itau = 1;
 222                 nwork = itau + *n;
 223
 224                 i__1 = *lwork - nwork + 1;
 225                 F77_FUNC(sgeqrf,SGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 226                         i__1, &ierr);
 227
 228                 i__1 = *n - 1;
 229                 i__2 = *n - 1;
 230                 F77_FUNC(slaset,SLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
 231                         lda);
 232                 ie = 1;
 233                 itauq = ie + *n;
 234                 itaup = itauq + *n;
 235                 nwork = itaup + *n;
 236
 237                 i__1 = *lwork - nwork + 1;
 238                 F77_FUNC(sgebrd,SGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
 239                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 240                 nwork = ie + *n;
 241
 242                 F77_FUNC(sbdsdc,SBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 243                          dum, idum, &work[nwork], &iwork[1], info);
 244
 245             } else {
 246                 iu = 1;
 247
 248                 ldwrku = *n;
 249                 itau = iu + ldwrku * *n;
 250                 nwork = itau + *n;
 251
 252                 i__1 = *lwork - nwork + 1;
 253                 F77_FUNC(sgeqrf,SGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 254                         i__1, &ierr);
 255                 F77_FUNC(slacpy,SLACPY)("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
 256
 257                 i__1 = *lwork - nwork + 1;
 258                 F77_FUNC(sorgqr,SORGQR)(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
 259                          &i__1, &ierr);
 260
 261                 i__1 = *n - 1;
 262                 i__2 = *n - 1;
 263                 F77_FUNC(slaset,SLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
 264                         lda);
 265                 ie = itau;
 266                 itauq = ie + *n;
 267                 itaup = itauq + *n;
 268                 nwork = itaup + *n;
 269
 270                 i__1 = *lwork - nwork + 1;
 271                 F77_FUNC(sgebrd,SGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
 272                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 273
 274                 F77_FUNC(sbdsdc,SBDSDC)("U", "I", n, &s[1], &work[ie], &work[iu], n, &vt[
 275                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 276                         info);
 277
 278                 i__1 = *lwork - nwork + 1;
 279                 F77_FUNC(sormbr,SORMBR)("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
 280                         itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
 281                         ierr);
 282                 i__1 = *lwork - nwork + 1;
 283                 F77_FUNC(sormbr,SORMBR)("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
 284                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
 285                         ierr);
 286
 287                 F77_FUNC(sgemm,SGEMM)("N", "N", m, n, n, &one, &u[u_offset], ldu, &work[iu]
 288                         , &ldwrku, &zero, &a[a_offset], lda);
 289
 290                 F77_FUNC(slacpy,SLACPY)("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
 291
 292             }
 293
 294         } else {
 295             ie = 1;
 296             itauq = ie + *n;
 297             itaup = itauq + *n;
 298             nwork = itaup + *n;
 299
 300             i__1 = *lwork - nwork + 1;
 301             F77_FUNC(sgebrd,SGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
 302                     work[itaup], &work[nwork], &i__1, &ierr);
 303             if (wntqn) {
 304
 305                 F77_FUNC(sbdsdc,SBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 306                          dum, idum, &work[nwork], &iwork[1], info);
 307             } else {
 308
 309                 F77_FUNC(slaset,SLASET)("F", m, m, &zero, &zero, &u[u_offset], ldu);
 310                 F77_FUNC(sbdsdc,SBDSDC)("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
 311                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 312                         info);
 313
 314                 i__1 = *m - *n;
 315                 i__2 = *m - *n;
 316                 F77_FUNC(slaset,SLASET)("F", &i__1, &i__2, &zero, &one, &u[*n + 1 + (*n +
 317                         1) * u_dim1], ldu);
 318
 319                 i__1 = *lwork - nwork + 1;
 320                 F77_FUNC(sormbr,SORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
 321                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 322                 i__1 = *lwork - nwork + 1;
 323                 F77_FUNC(sormbr,SORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
 324                         itaup], &vt[vt_offset],ldvt,&work[nwork],&i__1,&ierr);
 325             }
 326
 327         }
 328
 329     } else {
 330
 331         if (*n >= mnthr) {
 332
 333             if (wntqn) {
 334
 335                 itau = 1;
 336                 nwork = itau + *m;
 337
 338                 i__1 = *lwork - nwork + 1;
 339                 F77_FUNC(sgelqf,SGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 340                         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 = 1;
 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                 nwork = ie + *m;
 355
 356                 F77_FUNC(sbdsdc,SBDSDC)("U", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 357                          dum, idum, &work[nwork], &iwork[1], info);
 358
 359             } else {
 360
 361                 ivt = 1;
 362
 363                 ldwkvt = *m;
 364                 itau = ivt + ldwkvt * *m;
 365                 nwork = itau + *m;
 366
 367                 i__1 = *lwork - nwork + 1;
 368                 F77_FUNC(sgelqf,SGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 369                         i__1, &ierr);
 370                 F77_FUNC(slacpy,SLACPY)("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
 371
 372                 i__1 = *lwork - nwork + 1;
 373                 F77_FUNC(sorglq,SORGLQ)(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
 374                         nwork], &i__1, &ierr);
 375
 376                 i__1 = *m - 1;
 377                 i__2 = *m - 1;
 378                 F77_FUNC(slaset,SLASET)("U", &i__1, &i__2, &zero, &zero, &a[(a_dim1*2) +
 379                         1], lda);
 380                 ie = itau;
 381                 itauq = ie + *m;
 382                 itaup = itauq + *m;
 383                 nwork = itaup + *m;
 384
 385                 i__1 = *lwork - nwork + 1;
 386                 F77_FUNC(sgebrd,SGEBRD)(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
 387                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 388
 389                 F77_FUNC(sbdsdc,SBDSDC)("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
 390                         work[ivt], &ldwkvt, dum, idum, &work[nwork], &iwork[1]
 391                         , info);
 392
 393                 i__1 = *lwork - nwork + 1;
 394                 F77_FUNC(sormbr,SORMBR)("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
 395                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 396                 i__1 = *lwork - nwork + 1;
 397                 F77_FUNC(sormbr,SORMBR)("P", "R", "T", m, m, m, &a[a_offset], lda, &work[
 398                         itaup], &work[ivt], &ldwkvt, &work[nwork], &i__1, &
 399                         ierr);
 400
 401                 F77_FUNC(sgemm,SGEMM)("N", "N", m, n, m, &one, &work[ivt], &ldwkvt, &vt[
 402                         vt_offset], ldvt, &zero, &a[a_offset], lda);
 403
 404                 F77_FUNC(slacpy,SLACPY)("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
 405
 406             }
 407
 408         } else {
 409
 410             ie = 1;
 411             itauq = ie + *m;
 412             itaup = itauq + *m;
 413             nwork = itaup + *m;
 414
 415             i__1 = *lwork - nwork + 1;
 416             F77_FUNC(sgebrd,SGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
 417                     work[itaup], &work[nwork], &i__1, &ierr);
 418             if (wntqn) {
 419
 420                 F77_FUNC(sbdsdc,SBDSDC)("L", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 421                          dum, idum, &work[nwork], &iwork[1], info);
 422             } else {
 423                 F77_FUNC(slaset,SLASET)("F", n, n, &zero, &zero, &vt[vt_offset], ldvt);
 424                 F77_FUNC(sbdsdc,SBDSDC)("L", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &vt[
 425                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 426                         info);
 427
 428                 i__1 = *n - *m;
 429                 i__2 = *n - *m;
 430                 F77_FUNC(slaset,SLASET)("F", &i__1, &i__2, &zero, &one, &vt[*m + 1 + (*m +
 431                         1) * vt_dim1], ldvt);
 432
 433                 i__1 = *lwork - nwork + 1;
 434                 F77_FUNC(sormbr,SORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
 435                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 436                 i__1 = *lwork - nwork + 1;
 437                 F77_FUNC(sormbr,SORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
 438                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
 439                         ierr);
 440             }
 441
 442         }
 443
 444     }
 445
 446     if (iscl == 1) {
 447         if (anrm > bignum) {
 448             F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
 449                     minmn, &ierr);
 450         }
 451         if (anrm < smlnum) {
 452             F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
 453                     minmn, &ierr);
 454         }
 455     }
 456
 457     work[1] = (float) maxwrk;
 458
 459     return;
 460
 461 }
 462
 463