src/gromacs/linearalgebra/gmx_arpack.h

   1 /*
   2  *
   3  * This file is part of Gromacs        Copyright (c) 1991-2004
   4  * David van der Spoel, Erik Lindahl, University of Groningen.
   5  *
   6  * This file contains a subset of ARPACK functions to perform
   7  * diagonalization and SVD for sparse matrices in Gromacs.
   8  *
   9  * The code has been translated to C to avoid being dependent on
  10  * a Fotran compiler, and it has been made threadsafe by using
  11  * additional workspace arrays to store data during reverse communication.
  12  *
  13  * You might prefer the original ARPACK library for general use, but
  14  * in case you want to this version can be redistributed freely, just
  15  * as the original library. However, please make clear that it is the
  16  * hacked version from Gromacs so any bugs are blamed on us and not
  17  * the original authors. You should also be aware that the double
  18  * precision work array workd needs to be of size (3*N+4) here
  19  * (4 more than the general library), and there is an extra argument
  20  * iwork, which should be an integer work array of length 80.
  21  *
  22  * ARPACK was written by
  23  *
  24  *     Danny Sorensen               Phuong Vu
  25  *    Riconst chard Lehoucq              CRPC / Rice University
  26  *    Dept. of Computational &     Houston, Texas
  27  *    Applied Mathematics
  28  *    Rice University
  29  *    Houston, Texas
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Selected routines from ARPACK
  34  *
  35  * This file contains a subset of ARPACK functions to perform
  36  * diagonalization and SVD for sparse matrices in Gromacs.
  37  *
  38  * Consult the main ARPACK site for detailed documentation:
  39  * http://www.caam.rice.edu/software/ARPACK/
  40  *
  41  * Below, we just list the options and any specific differences
  42  * from ARPACK. The code is essentially the same, but the routines
  43  * have been made thread-safe by using extra workspace arrays.
  44  */
  45 #ifndef GMX_ARPACK_H
  46 #define GMX_ARPACK_H
  47
  48 #ifdef HAVE_CONFIG_H
  49 #include "config.h"
  50 #endif
  51 #ifndef F77_FUNC
  52 #define F77_FUNC(name,NAME) name ## _
  53 #endif
  54
  55 #ifdef __cplusplus
  56 extern "C" {
  57 #endif
  58
  59 /*! \brief Implicitly Restarted Arnoldi Iteration, double precision.
  60  *
  61  *  Reverse communication interface for the Implicitly Restarted Arnoldi
  62  *  Iteration.  For symmetric problems this reduces to a variant of the
  63  *  Lanczos method. See the ARPACK site for details.
  64  *
  65  *  \param ido     Reverse communication flag. Set to 0 first time.
  66  *                 Upon return with ido=-1 or ido=1 you should calculate
  67  *                 Y=A*X and recall the routine. Return with ido=2 means
  68  *                 Y=B*X should be calculated. ipntr[0] is the pointer in
  69  *                 workd for X, ipntr[1] is the index for Y.
  70  *                 Return with ido=99 means it finished.
  71  *  \param bmat    'I' for standard eigenproblem, 'G' for generalized.
  72  *  \param n       Order of eigenproblem.
  73  *  \param which   Which eigenvalues to calculate. 'LA' for largest
  74  *                 algebraic, 'SA' for smallest algebraic, 'LM' for largest
  75  *                 magnitude, 'SM' for smallest magnitude, and finally
  76  *                 'BE' (both ends) to calculate half from each end of
  77  *                 the spectrum.
  78  *  \param nev     Number of eigenvalues to calculate. 0<nev<n.
  79  *  \param tol     Tolerance. Machine precision of it is 0.
  80  *  \param resid   Optional starting residual vector at input if info=1,
  81  *                 otherwise a random one is used. Final residual vector on
  82  *                 return.
  83  *  \param ncv     Number of columns in matrix v.
  84  *  \param v       N*NCV matrix. V contain the Lanczos basis vectors.
  85  *  \param ldv     Leading dimension of v.
  86  *  \param iparam  Integer array, size 11. Same contents as arpack.
  87  *  \param ipntr   Integer array, size 11. Points to starting locations
  88  *                 in the workd/workl arrays. Same contents as arpack.
  89  *  \param workd   Double precision work array, length 3*n+4.
  90  *                 Provide the same array for all calls, and don't touch it.
  91  *                 IMPORTANT: This is 4 units larger than standard ARPACK!
  92  *  \param iwork   Integer work array, size 80.
  93  *                 Provide the same array for all calls, and don't touch it.
  94  *                 IMPORTANT: New argument compared to standard ARPACK!
  95  *  \param workl   Double precision work array, length lwork.
  96  *  \param lworkl  Length of the work array workl. Must be at least ncv*(ncv+8)
  97  *  \param info    Set info to 0 to use random initial residual vector,
  98  *                 or to 1 if you provide a one. On output, info=0 means
  99  *                 normal exit, 1 that max number of iterations was reached,
 100  *                 and 3 that no shifts could be applied. Negative numbers
 101  *                 correspond to errors in the arguments provided.
 102  */
 103 void
 104 F77_FUNC(dsaupd,DSAUPD)(int *     ido,
 105                         const char *    bmat,
 106                         int *     n,
 107                         const char *      which,
 108                         int *     nev,
 109                         double *  tol,
 110                         double *  resid,
 111                         int *     ncv,
 112                         double *  v,
 113                         int *     ldv,
 114                         int *     iparam,
 115                         int *     ipntr,
 116                         double *  workd,
 117                         int *     iwork,
 118                         double *  workl,
 119                         int *     lworkl,
 120                         int *     info);
 121
 122
 123
 124 /*! \brief Get eigenvalues/vectors after Arnoldi iteration, double prec.
 125  *
 126  *  See the ARPACK site for details. You must have finished the interative
 127  *  part with dsaupd() before calling this function.
 128  *
 129  *  \param rvec    1 if you want eigenvectors, 0 if not.
 130  *  \param howmny  'A' if you want all nvec vectors, 'S' if you
 131  *                 provide a subset selection in select[].
 132  *  \param select  Integer array, dimension nev. Indices of the
 133  *                 eigenvectors to calculate. Fortran code means we
 134  *                 start counting on 1. This array must be given even in
 135  *                 howmny is 'A'. (Arpack documentation is wrong on this).
 136  *  \param d       Double precision array, length nev. Eigenvalues.
 137  *  \param z       Double precision array, n*nev. Eigenvectors.
 138  *  \param ldz     Leading dimension of z. Normally n.
 139  *  \param sigma   Shift if iparam[6] is 3,4, or 5. Ignored otherwise.
 140  *  \param bmat    Provide the same argument as you did to dsaupd()
 141  *  \param n       Provide the same argument as you did to dsaupd()
 142  *  \param which   Provide the same argument as you did to dsaupd()
 143  *  \param nev     Provide the same argument as you did to dsaupd()
 144  *  \param tol     Provide the same argument as you did to dsaupd()
 145  *  \param resid   Provide the same argument as you did to dsaupd()
 146  *                 The array must not be touched between the two function calls!
 147  *  \param ncv     Provide the same argument as you did to dsaupd()
 148  *  \param v       Provide the same argument as you did to dsaupd()
 149  *                 The array must not be touched between the two function calls!
 150  *  \param ldv     Provide the same argument as you did to dsaupd()
 151  *  \param iparam  Provide the same argument as you did to dsaupd()
 152  *                 The array must not be touched between the two function calls!
 153  *  \param ipntr   Provide the same argument as you did to dsaupd()
 154  *                 The array must not be touched between the two function calls!
 155  *  \param workd   Provide the same argument as you did to dsaupd()
 156  *                 The array must not be touched between the two function calls!
 157  *  \param workl   Double precision work array, length lwork.
 158  *                 The array must not be touched between the two function calls!
 159  *  \param lworkl  Provide the same argument as you did to dsaupd()
 160  *  \param info    Provide the same argument as you did to dsaupd()
 161  */
 162 void
 163 F77_FUNC(dseupd,DSEUPD)(int *     rvec,
 164                         const char *    howmny,
 165                         int *     select,
 166                         double *  d,
 167                         double *  z,
 168                         int *     ldz,
 169                         double *  sigma,
 170                         const char *    bmat,
 171                         int *     n,
 172                         const char *    which,
 173                         int *     nev,
 174                         double *  tol,
 175                         double *  resid,
 176                         int *     ncv,
 177                         double *  v,
 178                         int *     ldv,
 179                         int *     iparam,
 180                         int *     ipntr,
 181                         double *  workd,
 182                         double *  workl,
 183                         int *     lworkl,
 184                         int *     info);
 185
 186
 187
 188
 189
 190 /*! \brief Implicitly Restarted Arnoldi Iteration, single precision.
 191  *
 192  *  Reverse communication interface for the Implicitly Restarted Arnoldi
 193  *  Iteration.  For symmetric problems this reduces to a variant of the
 194  *  Lanczos method. See the ARPACK site for details.
 195  *
 196  *  \param ido     Reverse communication flag. Set to 0 first time.
 197  *                 Upon return with ido=-1 or ido=1 you should calculate
 198  *                 Y=A*X and recall the routine. Return with ido=2 means
 199  *                 Y=B*X should be calculated. ipntr[0] is the pointer in
 200  *                 workd for X, ipntr[1] is the index for Y.
 201  *                 Return with ido=99 means it finished.
 202  *  \param bmat    'I' for standard eigenproblem, 'G' for generalized.
 203  *  \param n       Order of eigenproblem.
 204  *  \param which   Which eigenvalues to calculate. 'LA' for largest
 205  *                 algebraic, 'SA' for smallest algebraic, 'LM' for largest
 206  *                 magnitude, 'SM' for smallest magnitude, and finally
 207  *                 'BE' (both ends) to calculate half from each end of
 208  *                 the spectrum.
 209  *  \param nev     Number of eigenvalues to calculate. 0<nev<n.
 210  *  \param tol     Tolerance. Machine precision of it is 0.
 211  *  \param resid   Optional starting residual vector at input if info=1,
 212  *                 otherwise a random one is used. Final residual vector on
 213  *                 return.
 214  *  \param ncv     Number of columns in matrix v.
 215  *  \param v       N*NCV matrix. V contain the Lanczos basis vectors.
 216  *  \param ldv     Leading dimension of v.
 217  *  \param iparam  Integer array, size 11. Same contents as arpack.
 218  *  \param ipntr   Integer array, size 11. Points to starting locations
 219  *                 in the workd/workl arrays. Same contents as arpack.
 220  *  \param workd   Single precision work array, length 3*n+4.
 221  *                 Provide the same array for all calls, and don't touch it.
 222  *                 IMPORTANT: This is 4 units larger than standard ARPACK!
 223  *  \param iwork   Integer work array, size 80.
 224  *                 Provide the same array for all calls, and don't touch it.
 225  *                 IMPORTANT: New argument compared to standard ARPACK!
 226  *  \param workl   Single precision work array, length lwork.
 227  *  \param lworkl  Length of the work array workl. Must be at least ncv*(ncv+8)
 228  *  \param info    Set info to 0 to use random initial residual vector,
 229  *                 or to 1 if you provide a one. On output, info=0 means
 230  *                 normal exit, 1 that max number of iterations was reached,
 231  *                 and 3 that no shifts could be applied. Negative numbers
 232  *                 correspond to errors in the arguments provided.
 233  */
 234 void
 235 F77_FUNC(ssaupd,SSAUPD)(int *     ido,
 236                         const char *    bmat,
 237                         int *     n,
 238                         const char *    which,
 239                         int *     nev,
 240                         float *   tol,
 241                         float *   resid,
 242                         int *     ncv,
 243                         float *   v,
 244                         int *     ldv,
 245                         int *     iparam,
 246                         int *     ipntr,
 247                         float *   workd,
 248                         int *     iwork,
 249                         float *   workl,
 250                         int *     lworkl,
 251                         int *     info);
 252
 253
 254
 255
 256
 257 /*! \brief Get eigenvalues/vectors after Arnoldi iteration, single prec.
 258  *
 259  *  See the ARPACK site for details. You must have finished the interative
 260  *  part with ssaupd() before calling this function.
 261  *
 262  *  \param rvec    1 if you want eigenvectors, 0 if not.
 263  *  \param howmny  'A' if you want all nvec vectors, 'S' if you
 264  *                 provide a subset selection in select[].
 265  *  \param select  Integer array, dimension nev. Indices of the
 266  *                 eigenvectors to calculate. Fortran code means we
 267  *                 start counting on 1. This array must be given even in
 268  *                 howmny is 'A'. (Arpack documentation is wrong on this).
 269  *  \param d       Single precision array, length nev. Eigenvalues.
 270  *  \param z       Single precision array, n*nev. Eigenvectors.
 271  *  \param ldz     Leading dimension of z. Normally n.
 272  *  \param sigma   Shift if iparam[6] is 3,4, or 5. Ignored otherwise.
 273  *  \param bmat    Provide the same argument as you did to ssaupd()
 274  *  \param n       Provide the same argument as you did to ssaupd()
 275  *  \param which   Provide the same argument as you did to ssaupd()
 276  *  \param nev     Provide the same argument as you did to ssaupd()
 277  *  \param tol     Provide the same argument as you did to ssaupd()
 278  *  \param resid   Provide the same argument as you did to ssaupd()
 279  *                 The array must not be touched between the two function calls!
 280  *  \param ncv     Provide the same argument as you did to ssaupd()
 281  *  \param v       Provide the same argument as you did to ssaupd()
 282  *                 The array must not be touched between the two function calls!
 283  *  \param ldv     Provide the same argument as you did to ssaupd()
 284  *  \param iparam  Provide the same argument as you did to ssaupd()
 285  *                 The array must not be touched between the two function calls!
 286  *  \param ipntr   Provide the same argument as you did to ssaupd()
 287  *                 The array must not be touched between the two function calls!
 288  *  \param workd   Provide the same argument as you did to ssaupd()
 289  *                 The array must not be touched between the two function calls!
 290  *  \param workl   Single precision work array, length lwork.
 291  *                 The array must not be touched between the two function calls!
 292  *  \param lworkl  Provide the same argument as you did to ssaupd()
 293  *  \param info    Provide the same argument as you did to ssaupd()
 294  */
 295 void
 296 F77_FUNC(sseupd,SSEUPD)(int *     rvec,
 297                         const char *    howmny,
 298                         int *     select,
 299                         float *   d,
 300                         float *   z,
 301                         int *     ldz,
 302                         float *   sigma,
 303                         const char *    bmat,
 304                         int *     n,
 305                         const char *    which,
 306                         int *     nev,
 307                         float *   tol,
 308                         float *   resid,
 309                         int *     ncv,
 310                         float *   v,
 311                         int *     ldv,
 312                         int *     iparam,
 313                         int *     ipntr,
 314                         float *   workd,
 315                         float *   workl,
 316                         int *     lworkl,
 317                         int *     info);
 318
 319 #ifdef __cplusplus
 320 }
 321 #endif
 322
 323 #endif
 324