src/gromacs/linearalgebra/gmx_lapack/dgesdd.c

   1 #include <math.h>
   2 #include "types/simple.h"
   3
   4 #include "../gmx_blas.h"
   5 #include "../gmx_lapack.h"
   6 #include "lapack_limits.h"
   7
   8
   9 void
  10 F77_FUNC(dgesdd,DGESDD)(const char *jobz,
  11         int *m,
  12         int *n,
  13         double *a,
  14         int *lda,
  15         double *s,
  16         double *u,
  17         int *ldu,
  18         double *vt,
  19         int *ldvt,
  20         double *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     double dum[1], eps;
  29     int ivt, iscl;
  30     double 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     double bignum;
  38     int minwrk, ldwrku, maxwrk, ldwkvt;
  39     double smlnum,minval, safemin;
  40     int wntqas, lquery;
  41     int c__0 = 0;
  42     int c__1 = 1;
  43     double zero = 0.0;
  44     double 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] = (double) maxwrk;
 150     }
 151
 152
 153     if( lquery != 0)
 154     {
 155         return;
 156     }
 157
 158
 159     if (*m == 0 || *n == 0) {
 160         if (*lwork >= 1) {
 161             work[1] = 1.;
 162         }
 163         return;
 164     }
 165     eps = GMX_DOUBLE_EPS;
 166     minval = GMX_DOUBLE_MIN;
 167     safemin = minval / eps;
 168     smlnum = sqrt(safemin) / eps;
 169
 170
 171     bignum = 1. / smlnum;
 172
 173
 174     anrm = F77_FUNC(dlange,DLANGE)("M", m, n, &a[a_offset], lda, dum);
 175     iscl = 0;
 176     if (anrm > 0. && anrm < smlnum) {
 177         iscl = 1;
 178         F77_FUNC(dlascl,DLASCL)("G",&c__0,&c__0,&anrm,&smlnum,m,n,&a[a_offset],lda,&ierr);
 179     } else if (anrm > bignum) {
 180         iscl = 1;
 181         F77_FUNC(dlascl,DLASCL)("G",&c__0,&c__0,&anrm,&bignum,m,n,&a[a_offset],lda,&ierr);
 182     }
 183
 184     if (*m >= *n) {
 185         if (*m >= mnthr) {
 186
 187             if (wntqn) {
 188
 189                 itau = 1;
 190                 nwork = itau + *n;
 191
 192                 i__1 = *lwork - nwork + 1;
 193                 F77_FUNC(dgeqrf,DGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 194                         i__1, &ierr);
 195
 196                 i__1 = *n - 1;
 197                 i__2 = *n - 1;
 198                 F77_FUNC(dlaset,DLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
 199                         lda);
 200                 ie = 1;
 201                 itauq = ie + *n;
 202                 itaup = itauq + *n;
 203                 nwork = itaup + *n;
 204
 205                 i__1 = *lwork - nwork + 1;
 206                 F77_FUNC(dgebrd,DGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
 207                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 208                 nwork = ie + *n;
 209
 210                 F77_FUNC(dbdsdc,DBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 211                          dum, idum, &work[nwork], &iwork[1], info);
 212
 213             } else {
 214                 iu = 1;
 215
 216                 ldwrku = *n;
 217                 itau = iu + ldwrku * *n;
 218                 nwork = itau + *n;
 219
 220                 i__1 = *lwork - nwork + 1;
 221                 F77_FUNC(dgeqrf,DGEQRF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 222                         i__1, &ierr);
 223                 F77_FUNC(dlacpy,DLACPY)("L", m, n, &a[a_offset], lda, &u[u_offset], ldu);
 224
 225                 i__1 = *lwork - nwork + 1;
 226                 F77_FUNC(dorgqr,DORGQR)(m, m, n, &u[u_offset], ldu, &work[itau], &work[nwork],
 227                          &i__1, &ierr);
 228
 229                 i__1 = *n - 1;
 230                 i__2 = *n - 1;
 231                 F77_FUNC(dlaset,DLASET)("L", &i__1, &i__2, &zero, &zero, &a[a_dim1 + 2],
 232                         lda);
 233                 ie = itau;
 234                 itauq = ie + *n;
 235                 itaup = itauq + *n;
 236                 nwork = itaup + *n;
 237
 238                 i__1 = *lwork - nwork + 1;
 239                 F77_FUNC(dgebrd,DGEBRD)(n, n, &a[a_offset], lda, &s[1], &work[ie], &work[
 240                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 241
 242                 F77_FUNC(dbdsdc,DBDSDC)("U", "I", n, &s[1], &work[ie], &work[iu], n, &vt[
 243                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 244                         info);
 245
 246                 i__1 = *lwork - nwork + 1;
 247                 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", n, n, n, &a[a_offset], lda, &work[
 248                         itauq], &work[iu], &ldwrku, &work[nwork], &i__1, &
 249                         ierr);
 250                 i__1 = *lwork - nwork + 1;
 251                 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, n, &a[a_offset], lda, &work[
 252                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
 253                         ierr);
 254
 255                 F77_FUNC(dgemm,DGEMM)("N", "N", m, n, n, &one, &u[u_offset], ldu, &work[iu]
 256                         , &ldwrku, &zero, &a[a_offset], lda);
 257
 258                 F77_FUNC(dlacpy,DLACPY)("F", m, n, &a[a_offset], lda, &u[u_offset], ldu);
 259
 260             }
 261
 262         } else {
 263             ie = 1;
 264             itauq = ie + *n;
 265             itaup = itauq + *n;
 266             nwork = itaup + *n;
 267
 268             i__1 = *lwork - nwork + 1;
 269             F77_FUNC(dgebrd,DGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
 270                     work[itaup], &work[nwork], &i__1, &ierr);
 271             if (wntqn) {
 272
 273                 F77_FUNC(dbdsdc,DBDSDC)("U", "N", n, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 274                          dum, idum, &work[nwork], &iwork[1], info);
 275             } else {
 276
 277                 F77_FUNC(dlaset,DLASET)("F", m, m, &zero, &zero, &u[u_offset], ldu);
 278                 F77_FUNC(dbdsdc,DBDSDC)("U", "I", n, &s[1], &work[ie], &u[u_offset], ldu, &vt[
 279                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 280                         info);
 281
 282                 i__1 = *m - *n;
 283                 i__2 = *m - *n;
 284                 F77_FUNC(dlaset,DLASET)("F", &i__1, &i__2, &zero, &one, &u[*n + 1 + (*n +
 285                         1) * u_dim1], ldu);
 286
 287                 i__1 = *lwork - nwork + 1;
 288                 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
 289                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 290                 i__1 = *lwork - nwork + 1;
 291                 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
 292                         itaup], &vt[vt_offset],ldvt,&work[nwork],&i__1,&ierr);
 293             }
 294
 295         }
 296
 297     } else {
 298
 299         if (*n >= mnthr) {
 300
 301             if (wntqn) {
 302
 303                 itau = 1;
 304                 nwork = itau + *m;
 305
 306                 i__1 = *lwork - nwork + 1;
 307                 F77_FUNC(dgelqf,DGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 308                         i__1, &ierr);
 309
 310                 i__1 = *m - 1;
 311                 i__2 = *m - 1;
 312                 F77_FUNC(dlaset,DLASET)("U", &i__1, &i__2, &zero, &zero, &a[(a_dim1*2) +
 313                         1], lda);
 314                 ie = 1;
 315                 itauq = ie + *m;
 316                 itaup = itauq + *m;
 317                 nwork = itaup + *m;
 318
 319                 i__1 = *lwork - nwork + 1;
 320                 F77_FUNC(dgebrd,DGEBRD)(m, m, &a[a_offset], lda, &s[1], &work[ie], &work[
 321                         itauq], &work[itaup], &work[nwork], &i__1, &ierr);
 322                 nwork = ie + *m;
 323
 324                 F77_FUNC(dbdsdc,DBDSDC)("U", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 325                          dum, idum, &work[nwork], &iwork[1], info);
 326
 327             } else {
 328
 329                 ivt = 1;
 330
 331                 ldwkvt = *m;
 332                 itau = ivt + ldwkvt * *m;
 333                 nwork = itau + *m;
 334
 335                 i__1 = *lwork - nwork + 1;
 336                 F77_FUNC(dgelqf,DGELQF)(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &
 337                         i__1, &ierr);
 338                 F77_FUNC(dlacpy,DLACPY)("U", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
 339
 340                 i__1 = *lwork - nwork + 1;
 341                 F77_FUNC(dorglq,DORGLQ)(n, n, m, &vt[vt_offset], ldvt, &work[itau], &work[
 342                         nwork], &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 = itau;
 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
 357                 F77_FUNC(dbdsdc,DBDSDC)("U", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &
 358                         work[ivt], &ldwkvt, dum, idum, &work[nwork], &iwork[1]
 359                         , info);
 360
 361                 i__1 = *lwork - nwork + 1;
 362                 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, m, &a[a_offset], lda, &work[
 363                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 364                 i__1 = *lwork - nwork + 1;
 365                 F77_FUNC(dormbr,DORMBR)("P", "R", "T", m, m, m, &a[a_offset], lda, &work[
 366                         itaup], &work[ivt], &ldwkvt, &work[nwork], &i__1, &
 367                         ierr);
 368
 369                 F77_FUNC(dgemm,DGEMM)("N", "N", m, n, m, &one, &work[ivt], &ldwkvt, &vt[
 370                         vt_offset], ldvt, &zero, &a[a_offset], lda);
 371
 372                 F77_FUNC(dlacpy,DLACPY)("F", m, n, &a[a_offset], lda, &vt[vt_offset], ldvt);
 373
 374             }
 375
 376         } else {
 377
 378             ie = 1;
 379             itauq = ie + *m;
 380             itaup = itauq + *m;
 381             nwork = itaup + *m;
 382
 383             i__1 = *lwork - nwork + 1;
 384             F77_FUNC(dgebrd,DGEBRD)(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], &
 385                     work[itaup], &work[nwork], &i__1, &ierr);
 386             if (wntqn) {
 387
 388                 F77_FUNC(dbdsdc,DBDSDC)("L", "N", m, &s[1], &work[ie], dum, &c__1, dum, &c__1,
 389                          dum, idum, &work[nwork], &iwork[1], info);
 390             } else {
 391                 F77_FUNC(dlaset,DLASET)("F", n, n, &zero, &zero, &vt[vt_offset], ldvt);
 392                 F77_FUNC(dbdsdc,DBDSDC)("L", "I", m, &s[1], &work[ie], &u[u_offset], ldu, &vt[
 393                         vt_offset], ldvt, dum, idum, &work[nwork], &iwork[1],
 394                         info);
 395
 396                 i__1 = *n - *m;
 397                 i__2 = *n - *m;
 398                 F77_FUNC(dlaset,DLASET)("F", &i__1, &i__2, &zero, &one, &vt[*m + 1 + (*m +
 399                         1) * vt_dim1], ldvt);
 400
 401                 i__1 = *lwork - nwork + 1;
 402                 F77_FUNC(dormbr,DORMBR)("Q", "L", "N", m, m, n, &a[a_offset], lda, &work[
 403                         itauq], &u[u_offset], ldu, &work[nwork], &i__1, &ierr);
 404                 i__1 = *lwork - nwork + 1;
 405                 F77_FUNC(dormbr,DORMBR)("P", "R", "T", n, n, m, &a[a_offset], lda, &work[
 406                         itaup], &vt[vt_offset], ldvt, &work[nwork], &i__1, &
 407                         ierr);
 408             }
 409
 410         }
 411
 412     }
 413
 414     if (iscl == 1) {
 415         if (anrm > bignum) {
 416             F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], &
 417                     minmn, &ierr);
 418         }
 419         if (anrm < smlnum) {
 420             F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], &
 421                     minmn, &ierr);
 422         }
 423     }
 424
 425     work[1] = (double) maxwrk;
 426
 427     return;
 428
 429 }
 430
 431