src/gmxlib/gmx_lapack/dgesdd.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013, 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 #include "gmx_blas.h"
  39 #include "gmx_lapack.h"
  40 #include "lapack_limits.h"
  41
  42
  43 void
  44 F77_FUNC(dgesdd,DGESDD)(const char *jobz,
  45         int *m,
  46         int *n,
  47         double *a,
  48         int *lda,
  49         double *s,
  50         double *u,
  51         int *ldu,
  52         double *vt,
  53         int *ldvt,
  54         double *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     double dum[1], eps;
  63     int ivt, iscl;
  64     double 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     double bignum;
  72     int minwrk, ldwrku, maxwrk, ldwkvt;
  73     double smlnum,minval, safemin;
  74     int wntqas, lquery;
  75     int c__0 = 0;
  76     int c__1 = 1;
  77     double zero = 0.0;
  78     double 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] = (double) maxwrk;
 184     }
 185
 186
 187     if( lquery != 0)
 188     {
 189         return;
 190     }
 191
 192
 193     if (*m == 0 || *n == 0) {
 194         if (*lwork >= 1) {
 195             work[1] = 1.;
 196         }
 197         return;
 198     }
 199     eps = GMX_DOUBLE_EPS;
 200     minval = GMX_DOUBLE_MIN;
 201     safemin = minval / eps;
 202     smlnum = sqrt(safemin) / eps;
 203
 204
 205     bignum = 1. / smlnum;
 206
 207
 208     anrm = F77_FUNC(dlange,DLANGE)("M", m, n, &a[a_offset], lda, dum);
 209     iscl = 0;
 210     if (anrm > 0. && anrm < smlnum) {
 211         iscl = 1;
 212         F77_FUNC(dlascl,DLASCL)("G",&c__0,&c__0,&anrm,&smlnum,m,n,&a[a_offset],lda,&ierr);
 213     } else if (anrm > bignum) {
 214         iscl = 1;
 215         F77_FUNC(dlascl,DLASCL)("G",&c__0,&c__0,&anrm,&bignum,m,n,&a[a_offset],lda,&ierr);
 216     }
 217
 218     if (*m >= *n) {
 219         if (*m >= mnthr) {
 220
 221             if (wntqn) {
 222
 223                 itau = 1;
 224                 nwork = itau + *n;
 225
 226                 i__1 = *lwork - nwork + 1;
 227                 F77_FUNC(dgeqrf,DGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 228                         i__1, &ierr);
 229
 230                 i__1 = *n - 1;
 231                 i__2 = *n - 1;
 232                 F77_FUNC(dlaset,DLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
 233                         lda);
 234                 ie = 1;
 235                 itauq = ie + *n;
 236                 itaup = itauq + *n;
 237                 nwork = itaup + *n;
 238
 239                 i__1 = *lwork - nwork + 1;
 240                 F77_FUNC(dgebrd,DGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
 241                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 242                 nwork = ie + *n;
 243
 244                 F77_FUNC(dbdsdc,DBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 245                          dum, idum, &work[nwork], &iwork[1], info);
 246
 247             } else {
 248                 iu = 1;
 249
 250                 ldwrku = *n;
 251                 itau = iu + ldwrku * *n;
 252                 nwork = itau + *n;
 253
 254                 i__1 = *lwork - nwork + 1;
 255                 F77_FUNC(dgeqrf,DGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 256                         i__1, &ierr);
 257                 F77_FUNC(dlacpy,DLACPY)("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
 258
 259                 i__1 = *lwork - nwork + 1;
 260                 F77_FUNC(dorgqr,DORGQR)(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
 261                          &i__1, &ierr);
 262
 263                 i__1 = *n - 1;
 264                 i__2 = *n - 1;
 265                 F77_FUNC(dlaset,DLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
 266                         lda);
 267                 ie = itau;
 268                 itauq = ie + *n;
 269                 itaup = itauq + *n;
 270                 nwork = itaup + *n;
 271
 272                 i__1 = *lwork - nwork + 1;
 273                 F77_FUNC(dgebrd,DGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
 274                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 275
 276                 F77_FUNC(dbdsdc,DBDSDC)("U", "I", n, &s[1], &work[ie], &work[iu], n, &vt[
 277                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 278                         info);
 279
 280                 i__1 = *lwork - nwork + 1;
 281                 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
 282                         itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
 283                         ierr);
 284                 i__1 = *lwork - nwork + 1;
 285                 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
 286                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
 287                         ierr);
 288
 289                 F77_FUNC(dgemm,DGEMM)("N", "N", m, n, n, &one, &u[u_offset], ldu, &work[iu]
 290                         , &ldwrku, &zero, &a[a_offset], lda);
 291
 292                 F77_FUNC(dlacpy,DLACPY)("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
 293
 294             }
 295
 296         } else {
 297             ie = 1;
 298             itauq = ie + *n;
 299             itaup = itauq + *n;
 300             nwork = itaup + *n;
 301
 302             i__1 = *lwork - nwork + 1;
 303             F77_FUNC(dgebrd,DGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
 304                     work[itaup], &work[nwork], &i__1, &ierr);
 305             if (wntqn) {
 306
 307                 F77_FUNC(dbdsdc,DBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 308                          dum, idum, &work[nwork], &iwork[1], info);
 309             } else {
 310
 311                 F77_FUNC(dlaset,DLASET)("F", m, m, &zero, &zero, &u[u_offset], ldu);
 312                 F77_FUNC(dbdsdc,DBDSDC)("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
 313                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 314                         info);
 315
 316                 i__1 = *m - *n;
 317                 i__2 = *m - *n;
 318                 F77_FUNC(dlaset,DLASET)("F", &i__1, &i__2, &zero, &one, &u[*n + 1 + (*n +
 319                         1) * u_dim1], ldu);
 320
 321                 i__1 = *lwork - nwork + 1;
 322                 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
 323                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 324                 i__1 = *lwork - nwork + 1;
 325                 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
 326                         itaup], &vt[vt_offset],ldvt,&work[nwork],&i__1,&ierr);
 327             }
 328
 329         }
 330
 331     } else {
 332
 333         if (*n >= mnthr) {
 334
 335             if (wntqn) {
 336
 337                 itau = 1;
 338                 nwork = itau + *m;
 339
 340                 i__1 = *lwork - nwork + 1;
 341                 F77_FUNC(dgelqf,DGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 342                         i__1, &ierr);
 343
 344                 i__1 = *m - 1;
 345                 i__2 = *m - 1;
 346                 F77_FUNC(dlaset,DLASET)("U", &i__1, &i__2, &zero, &zero, &a[(a_dim1*2) +
 347                         1], lda);
 348                 ie = 1;
 349                 itauq = ie + *m;
 350                 itaup = itauq + *m;
 351                 nwork = itaup + *m;
 352
 353                 i__1 = *lwork - nwork + 1;
 354                 F77_FUNC(dgebrd,DGEBRD)(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
 355                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 356                 nwork = ie + *m;
 357
 358                 F77_FUNC(dbdsdc,DBDSDC)("U", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 359                          dum, idum, &work[nwork], &iwork[1], info);
 360
 361             } else {
 362
 363                 ivt = 1;
 364
 365                 ldwkvt = *m;
 366                 itau = ivt + ldwkvt * *m;
 367                 nwork = itau + *m;
 368
 369                 i__1 = *lwork - nwork + 1;
 370                 F77_FUNC(dgelqf,DGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 371                         i__1, &ierr);
 372                 F77_FUNC(dlacpy,DLACPY)("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
 373
 374                 i__1 = *lwork - nwork + 1;
 375                 F77_FUNC(dorglq,DORGLQ)(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
 376                         nwork], &i__1, &ierr);
 377
 378                 i__1 = *m - 1;
 379                 i__2 = *m - 1;
 380                 F77_FUNC(dlaset,DLASET)("U", &i__1, &i__2, &zero, &zero, &a[(a_dim1*2) +
 381                         1], lda);
 382                 ie = itau;
 383                 itauq = ie + *m;
 384                 itaup = itauq + *m;
 385                 nwork = itaup + *m;
 386
 387                 i__1 = *lwork - nwork + 1;
 388                 F77_FUNC(dgebrd,DGEBRD)(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
 389                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 390
 391                 F77_FUNC(dbdsdc,DBDSDC)("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
 392                         work[ivt], &ldwkvt, dum, idum, &work[nwork], &iwork[1]
 393                         , info);
 394
 395                 i__1 = *lwork - nwork + 1;
 396                 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
 397                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 398                 i__1 = *lwork - nwork + 1;
 399                 F77_FUNC(dormbr,DORMBR)("P", "R", "T", m, m, m, &a[a_offset], lda, &work[
 400                         itaup], &work[ivt], &ldwkvt, &work[nwork], &i__1, &
 401                         ierr);
 402
 403                 F77_FUNC(dgemm,DGEMM)("N", "N", m, n, m, &one, &work[ivt], &ldwkvt, &vt[
 404                         vt_offset], ldvt, &zero, &a[a_offset], lda);
 405
 406                 F77_FUNC(dlacpy,DLACPY)("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
 407
 408             }
 409
 410         } else {
 411
 412             ie = 1;
 413             itauq = ie + *m;
 414             itaup = itauq + *m;
 415             nwork = itaup + *m;
 416
 417             i__1 = *lwork - nwork + 1;
 418             F77_FUNC(dgebrd,DGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
 419                     work[itaup], &work[nwork], &i__1, &ierr);
 420             if (wntqn) {
 421
 422                 F77_FUNC(dbdsdc,DBDSDC)("L", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 423                          dum, idum, &work[nwork], &iwork[1], info);
 424             } else {
 425                 F77_FUNC(dlaset,DLASET)("F", n, n, &zero, &zero, &vt[vt_offset], ldvt);
 426                 F77_FUNC(dbdsdc,DBDSDC)("L", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &vt[
 427                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 428                         info);
 429
 430                 i__1 = *n - *m;
 431                 i__2 = *n - *m;
 432                 F77_FUNC(dlaset,DLASET)("F", &i__1, &i__2, &zero, &one, &vt[*m + 1 + (*m +
 433                         1) * vt_dim1], ldvt);
 434
 435                 i__1 = *lwork - nwork + 1;
 436                 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
 437                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 438                 i__1 = *lwork - nwork + 1;
 439                 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
 440                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
 441                         ierr);
 442             }
 443
 444         }
 445
 446     }
 447
 448     if (iscl == 1) {
 449         if (anrm > bignum) {
 450             F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
 451                     minmn, &ierr);
 452         }
 453         if (anrm < smlnum) {
 454             F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
 455                     minmn, &ierr);
 456         }
 457     }
 458
 459     work[1] = (double) maxwrk;
 460
 461     return;
 462
 463 }
 464
 465