EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
misc_la.c
Go to the documentation of this file.
1 //============================================================================
2 // Routines for computing eigenvalues of a symmetric tridiagonal matrix.
3 // They are wrappers of the LAPACK routine DSTEV() or sstev_()
4 //============================================================================
5 
6 #include <stdio.h>
7 #include <string.h> // for memcpy, strcmp, strncmp, etc.
8 #include "def.h"
9 #include "blaslapack.h"
10 #include "struct.h"
11 #include "internal_proto.h"
12 
13 /* in Gram-Schmidt, use BLAS-2 GEMV or BLAS-1 DOT/AXPY */
14 #define USE_DGEMV 1
15 
36 int SymmTridEig(double *eigVal, double *eigVec, int n,
37  const double *diag, const double *sdiag) {
38  double tms = evsl_timer();
39  // compute eigenvalues and eigenvectors or eigvalues only
40  char jobz = eigVec ? 'V' : 'N';
41  int nn = n;
42  int ldz = n;
43  int info; // output flag
44  // copy diagonal and subdiagonal elements to alp and bet
45  double *alp = eigVal;
46  double *bet;
47  Malloc(bet, n-1, double);
48  memcpy(alp, diag, n*sizeof(double));
49  memcpy(bet, sdiag, (n-1)*sizeof(double));
50  // allocate storage for computation
51  double *sv = eigVec;
52  double *work = NULL;
53  if (jobz == 'V') {
54  Malloc(work, 2*n-2, double);
55  }
56  DSTEV(&jobz, &nn, alp, bet, sv, &ldz, work, &info);
57  // free memory
58  free(bet);
59  if (work) {
60  free(work);
61  }
62  if (info) {
63  printf("DSTEV ERROR: INFO %d\n", info);
64  save_vec(n, diag, "alp");
65  save_vec(n-1, sdiag, "bet");
66  exit(0);
67  }
68  double tme = evsl_timer();
69  evslstat.t_eig += tme - tms;
70  // return info
71  return info;
72 }
73 
93 int SymmTridEigS(double *eigVal, double *eigVec, int n, double vl, double vu,
94  int *nevO, const double *diag, const double *sdiag) {
95  double tms = evsl_timer();
96  char jobz = 'V'; // compute eigenvalues and eigenvectors
97  char range = 'V'; // compute eigenvalues in an interval
98 
99  // this does not use mwlapack for mex files
100  int info;
101  //int idum = 0;
102  //-------------------- isuppz not needed
103  int *isuppz;
104  Malloc(isuppz, 2*n, int);
105  //-------------------- real work array
106  double *work;
107  int lwork = 18*n;
108  Malloc(work, lwork, double);
109  //-------------------- int work array
110  int *iwork;
111  int liwork = 10*n;
112  Calloc(iwork, liwork, int);
113  //-------------------- copy diagonal + subdiagonal elements
114  // to alp and bet
115  double *alp;
116  double *bet;
117  Malloc(bet, n, double);
118  Malloc(alp, n, double);
119  //
120  memcpy(alp, diag, n*sizeof(double));
121  if (n > 1) {
122  memcpy(bet, sdiag, (n-1)*sizeof(double));
123  }
124 
125  //-------------------- allocate storage for computation
126  //logical tryrac = 1;
127  int tryrac = 1;
128  double t0 = vl;
129  double t1 = vu;
130 
131  DSTEMR(&jobz, &range, &n, alp, bet, &t0, &t1, NULL, NULL, nevO,
132  eigVal, eigVec, &n, &n, isuppz, &tryrac, work, &lwork,
133  iwork, &liwork, &info);
134 
135  if (info) {
136  printf("dstemr_ error %d\n", info);
137  }
138 
139  //-------------------- free memory
140  free(bet);
141  free(alp);
142  free(work);
143  free(iwork);
144  free(isuppz);
145 
146  double tme = evsl_timer();
147  evslstat.t_eig += tme - tms;
148  //
149  return info;
150 }
151 
152 
156 void SymEigenSolver(int n, double *A, int lda, double *Q, int ldq, double *lam) {
157  double tms = evsl_timer();
158  /* compute eigenvalues/vectors of A that n x n, symmetric
159  * eigenvalues saved in lam: the eigenvalues in ascending order
160  * eigenvectors saved in Q */
161  char jobz='V';/* want eigenvectors */
162  char uplo='U';/* use upper triangular part of the matrix */
163  /* copy A to Q */
164  int i,j;
165  for (i=0; i<n; i++) {
166  for (j=0; j<n; j++) {
167  Q[j+i*ldq] = A[j+i*lda];
168  }
169  }
170  /* first call: workspace query */
171  double work_size;
172  int lwork = -1, info;
173  DSYEV(&jobz, &uplo, &n, Q, &ldq, lam, &work_size, &lwork, &info);
174  if (info != 0) {
175  fprintf(stdout, "DSYEV error [query call]: %d\n", info);
176  exit(0);
177  }
178  /* second call: do the computation */
179  lwork = (int) work_size;
180  double *work;
181  Malloc(work, lwork, double);
182  DSYEV(&jobz, &uplo, &n, Q, &ldq, lam, work, &lwork, &info);
183  if (info != 0) {
184  fprintf(stdout, "DSYEV error [comp call]: %d\n", info);
185  exit(0);
186  }
187  free(work);
188  double tme = evsl_timer();
189  evslstat.t_eig += tme - tms;
190 }
191 
195 void CGS_DGKS(int n, int k, int i_max, double *Q, double *v, double *nrmv, double *w) {
196  double tms = evsl_timer();
197  double eta = 1.0 / sqrt(2.0);
198  int i, one=1;
199 #if USE_DGEMV
200  char cT = 'T', cN = 'N';
201  double done=1.0, dmone=-1.0, dzero=0.0;
202 #endif
203  double old_nrm = DNRM2(&n, v, &one);
204  double new_nrm = 0.0;
205 
206  for (i=0; i<i_max; i++) {
207 #if USE_DGEMV
208  DGEMV(&cT, &n, &k, &done, Q, &n, v, &one, &dzero, w, &one);
209  DGEMV(&cN, &n, &k, &dmone, Q, &n, w, &one, &done, v, &one);
210 #else
211  int j;
212  for (j=0; j<k; j++) {
213  double t = -DDOT(&n, &Q[j*n], &one, v, &one);
214  DAXPY(&n, &t, &Q[j*n], &one, v, &one);
215  }
216 #endif
217  new_nrm = DNRM2(&n, v, &one);
218  if (new_nrm > eta * old_nrm) {
219  break;
220  }
221  old_nrm = new_nrm;
222  }
223 
224  if (nrmv) {
225  *nrmv = new_nrm;
226  }
227  double tme = evsl_timer();
228  evslstat.t_reorth += tme - tms;
229 }
230 
235 void CGS_DGKS2(int n, int k, int i_max, double *Z, double *Q,
236  double *v, double *w) {
237  double tms = evsl_timer();
238  int i, one=1;
239 #if USE_DGEMV
240  char cT = 'T', cN = 'N';
241  double done=1.0, dmone=-1.0, dzero=0.0;
242 #endif
243  for (i=0; i<i_max; i++) {
244 #if USE_DGEMV
245  DGEMV(&cT, &n, &k, &done, Q, &n, v, &one, &dzero, w, &one);
246  DGEMV(&cN, &n, &k, &dmone, Z, &n, w, &one, &done, v, &one);
247 #else
248  int j;
249  for (j=0; j<k; j++) {
250  double t = -DDOT(&n, &Q[j*n], &one, v, &one);
251  DAXPY(&n, &t, &Z[j*n], &one, v, &one);
252  }
253 #endif
254  }
255  double tme = evsl_timer();
256  evslstat.t_reorth += tme - tms;
257 }
258 
259 // max number of reorthogonalizations
260 #define NGS_MAX 2
261 
269 void orth(double *V, int n, int k, double *Vo, double *work) {
270  int i;
271  int one=1;
272  int nk = n*k;
273  DCOPY(&nk, V, &one, Vo, &one);
274  double tt = DDOT(&n, Vo, &one, Vo, &one);
275  double nrmv = sqrt(tt);
276  double t = 1.0 / nrmv;
277  DSCAL(&n, &t, Vo, &one);
278  for (i = 1; i < k; i++) {
279  int istart = i*n;
280  CGS_DGKS(n, i, NGS_MAX, Vo, Vo+istart, &nrmv, work);
281  t = 1.0 / nrmv;
282  DSCAL(&n, &t, Vo+istart, &one);
283  }
284 }
285 
void SymEigenSolver(int n, double *A, int lda, double *Q, int ldq, double *lam)
interface to LAPACK SYMMETRIC EIGEN-SOLVER
Definition: misc_la.c:156
void DSTEV(char *jobz, int *n, double *diagonal, double *subdiagonal, double *V, int *ldz, double *work, int *info)
defs in EVSL
double DDOT(int *n, double *x, int *incx, double *y, int *incy)
void CGS_DGKS(int n, int k, int i_max, double *Q, double *v, double *nrmv, double *w)
Classical GS reortho with Daniel, Gragg, Kaufman, Stewart test.
Definition: misc_la.c:195
void save_vec(int n, const double *x, const char fn[])
Definition: dumps.c:38
void DSYEV(char *jobz, char *uplo, int *n, double *fa, int *lda, double *w, double *work, int *lwork, int *info)
void DSTEMR(char *jobz, char *range, int *n, double *D, double *E, double *VL, double *VU, int *IL, int *IU, int *M, double *W, double *Z, int *LDZ, int *NZC, int *ISUPPZ, logical *TRYRAC, double *WORK, int *LWORK, int *IWORK, int *LIWORK, int *INFO)
#define NGS_MAX
Definition: misc_la.c:260
void DSCAL(int *n, double *a, double *x, int *incx)
#define Calloc(base, nmem, type)
Definition: def.h:32
int SymmTridEig(double *eigVal, double *eigVec, int n, const double *diag, const double *sdiag)
compute all eigenvalues and eigenvectors of a symmetric tridiagonal matrix
Definition: misc_la.c:36
double t_reorth
Definition: struct.h:207
This file contains function prototypes and constant definitions internally used in EVSL...
void orth(double *V, int n, int k, double *Vo, double *work)
Orthogonalize columns of n-by-k matrix V.
Definition: misc_la.c:269
#define Malloc(base, nmem, type)
Definition: def.h:22
void DGEMV(char *trans, int *m, int *n, double *alpha, double *a, int *lda, double *x, int *incx, double *beta, double *y, int *incy)
evslStat evslstat
global statistics of EVSL
Definition: evsl.c:21
void CGS_DGKS2(int n, int k, int i_max, double *Z, double *Q, double *v, double *w)
Classical GS reortho. No test. just do i_max times used in generalized ev problems.
Definition: misc_la.c:235
structs used in evsl
void DAXPY(int *n, double *alpha, double *x, int *incx, double *y, int *incy)
double DNRM2(int *n, double *x, int *incx)
Defs for blaslapack routines.
int SymmTridEigS(double *eigVal, double *eigVec, int n, double vl, double vu, int *nevO, const double *diag, const double *sdiag)
compute eigenvalues and eigenvectors of a symmetric tridiagonal matrix in a slice ...
Definition: misc_la.c:93
double t_eig
Definition: struct.h:208
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
void DCOPY(int *n, double *dx, int *incx, double *dy, int *incy)