MAGMA  magma-1.4.0
Matrix Algebra on GPU and Multicore Architectures
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
ssyevd.cpp File Reference
#include "common_magma.h"
Include dependency graph for ssyevd.cpp:

Go to the source code of this file.

Functions

magma_int_t magma_ssyevd (char jobz, char uplo, magma_int_t n, float *a, magma_int_t lda, float *w, float *work, magma_int_t lwork, magma_int_t *iwork, magma_int_t liwork, magma_int_t *info)
 

Function Documentation

magma_int_t magma_ssyevd ( char  jobz,
char  uplo,
magma_int_t  n,
float *  a,
magma_int_t  lda,
float *  w,
float *  work,
magma_int_t  lwork,
magma_int_t iwork,
magma_int_t  liwork,
magma_int_t info 
)

Definition at line 17 of file ssyevd.cpp.

References __func__, blasf77_sscal(), dwork, get_current_time(), GetTimerValue(), lapackf77_lsame, lapackf77_slacpy(), lapackf77_slamch, lapackf77_slansy(), lapackf77_slascl(), lapackf77_ssterf, lapackf77_ssyevd(), MAGMA_ERR_DEVICE_ALLOC, magma_free, magma_get_ssytrd_nb(), magma_smalloc(), magma_sormtr(), magma_ssqrt, magma_sstedx(), magma_ssytrd(), MAGMA_SUCCESS, magma_xerbla(), MagmaLeft, MagmaLowerStr, MagmaNoTrans, MagmaNoVecStr, MagmaUpperStr, MagmaVecStr, and max.

24 {
25 /* -- MAGMA (version 1.4.0) --
26  Univ. of Tennessee, Knoxville
27  Univ. of California, Berkeley
28  Univ. of Colorado, Denver
29  August 2013
30 
31  Purpose
32  =======
33  SSYEVD computes all eigenvalues and, optionally, eigenvectors of
34  a real symmetric matrix A. If eigenvectors are desired, it uses a
35  divide and conquer algorithm.
36 
37  The divide and conquer algorithm makes very mild assumptions about
38  floating point arithmetic. It will work on machines with a guard
39  digit in add/subtract, or on those binary machines without guard
40  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
41  Cray-2. It could conceivably fail on hexadecimal or decimal machines
42  without guard digits, but we know of none.
43 
44  Arguments
45  =========
46  JOBZ (input) CHARACTER*1
47  = 'N': Compute eigenvalues only;
48  = 'V': Compute eigenvalues and eigenvectors.
49 
50  UPLO (input) CHARACTER*1
51  = 'U': Upper triangle of A is stored;
52  = 'L': Lower triangle of A is stored.
53 
54  N (input) INTEGER
55  The order of the matrix A. N >= 0.
56 
57  A (input/output) REAL array, dimension (LDA, N)
58  On entry, the symmetric matrix A. If UPLO = 'U', the
59  leading N-by-N upper triangular part of A contains the
60  upper triangular part of the matrix A. If UPLO = 'L',
61  the leading N-by-N lower triangular part of A contains
62  the lower triangular part of the matrix A.
63  On exit, if JOBZ = 'V', then if INFO = 0, A contains the
64  orthonormal eigenvectors of the matrix A.
65  If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
66  or the upper triangle (if UPLO='U') of A, including the
67  diagonal, is destroyed.
68 
69  LDA (input) INTEGER
70  The leading dimension of the array A. LDA >= max(1,N).
71 
72  W (output) DOUBLE PRECISION array, dimension (N)
73  If INFO = 0, the eigenvalues in ascending order.
74 
75  WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
76  On exit, if INFO = 0, WORK[0] returns the optimal LWORK.
77 
78  LWORK (input) INTEGER
79  The length of the array WORK.
80  If N <= 1, LWORK >= 1.
81  If JOBZ = 'N' and N > 1, LWORK >= 2*N + N*NB.
82  If JOBZ = 'V' and N > 1, LWORK >= max( 2*N + N*NB, 1 + 6*N + 2*N**2 ).
83  NB can be obtained through magma_get_ssytrd_nb(N).
84 
85  If LWORK = -1, then a workspace query is assumed; the routine
86  only calculates the optimal sizes of the WORK and IWORK
87  arrays, returns these values as the first entries of the WORK
88  and IWORK arrays, and no error message related to LWORK or
89  LIWORK is issued by XERBLA.
90 
91  IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
92  On exit, if INFO = 0, IWORK[0] returns the optimal LIWORK.
93 
94  LIWORK (input) INTEGER
95  The dimension of the array IWORK.
96  If N <= 1, LIWORK >= 1.
97  If JOBZ = 'N' and N > 1, LIWORK >= 1.
98  If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N.
99 
100  If LIWORK = -1, then a workspace query is assumed; the
101  routine only calculates the optimal sizes of the WORK and
102  IWORK arrays, returns these values as the first entries of
103  the WORK and IWORK arrays, and no error message related to
104  LWORK or LIWORK is issued by XERBLA.
105 
106  INFO (output) INTEGER
107  = 0: successful exit
108  < 0: if INFO = -i, the i-th argument had an illegal value
109  > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
110  to converge; i off-diagonal elements of an intermediate
111  tridiagonal form did not converge to zero;
112  if INFO = i and JOBZ = 'V', then the algorithm failed
113  to compute an eigenvalue while working on the submatrix
114  lying in rows and columns INFO/(N+1) through
115  mod(INFO,N+1).
116 
117  Further Details
118  ===============
119  Based on contributions by
120  Jeff Rutter, Computer Science Division, University of California
121  at Berkeley, USA
122 
123  Modified description of INFO. Sven, 16 Feb 05.
124  ===================================================================== */
125 
126  char uplo_[2] = {uplo, 0};
127  char jobz_[2] = {jobz, 0};
128  magma_int_t ione = 1;
129  magma_int_t izero = 0;
130  float d_one = 1.;
131 
132  float d__1;
133 
134  float eps;
135  magma_int_t inde;
136  float anrm;
137  float rmin, rmax;
138  float sigma;
139  magma_int_t iinfo, lwmin;
140  magma_int_t lower;
141  magma_int_t wantz;
142  magma_int_t indwk2, llwrk2;
143  magma_int_t iscale;
144  float safmin;
145  float bignum;
146  magma_int_t indtau;
147  magma_int_t indwrk, liwmin;
148  magma_int_t llwork;
149  float smlnum;
150  magma_int_t lquery;
151 
152  float* dwork;
153 
154  wantz = lapackf77_lsame(jobz_, MagmaVecStr);
155  lower = lapackf77_lsame(uplo_, MagmaLowerStr);
156  lquery = lwork == -1 || liwork == -1;
157 
158  *info = 0;
159 
160  if (! (wantz || lapackf77_lsame(jobz_, MagmaNoVecStr))) {
161  *info = -1;
162  } else if (! (lower || lapackf77_lsame(uplo_, MagmaUpperStr))) {
163  *info = -2;
164  } else if (n < 0) {
165  *info = -3;
166  } else if (lda < max(1,n)) {
167  *info = -5;
168  }
169 
171  if ( n <= 1 ) {
172  lwmin = 1;
173  liwmin = 1;
174  }
175  else if ( wantz ) {
176  lwmin = max( 2*n + n*nb, 1 + 6*n + 2*n*n );
177  liwmin = 3 + 5*n;
178  }
179  else {
180  lwmin = 2*n + n*nb;
181  liwmin = 1;
182  }
183  // multiply by 1+eps to ensure length gets rounded up,
184  // if it cannot be exactly represented in floating point.
185  work[0] = lwmin * (1. + lapackf77_slamch("Epsilon"));
186  iwork[0] = liwmin;
187 
188  if ((lwork < lwmin) && !lquery) {
189  *info = -8;
190  } else if ((liwork < liwmin) && ! lquery) {
191  *info = -10;
192  }
193 
194  if (*info != 0) {
195  magma_xerbla( __func__, -(*info) );
196  return *info;
197  }
198  else if (lquery) {
199  return *info;
200  }
201 
202  /* Quick return if possible */
203  if (n == 0) {
204  return *info;
205  }
206 
207  if (n == 1) {
208  w[0] = a[0];
209  if (wantz) {
210  a[0] = 1.;
211  }
212  return *info;
213  }
214 
215  /* Check if matrix is very small then just call LAPACK on CPU, no need for GPU */
216  if (n <= 128){
217  #ifdef ENABLE_DEBUG
218  printf("--------------------------------------------------------------\n");
219  printf(" warning matrix too small N=%d NB=%d, calling lapack on CPU \n", (int) n, (int) nb);
220  printf("--------------------------------------------------------------\n");
221  #endif
222  lapackf77_ssyevd(jobz_, uplo_,
223  &n, a, &lda,
224  w, work, &lwork,
225  iwork, &liwork, info);
226  return *info;
227  }
228 
229  /* Get machine constants. */
230  safmin = lapackf77_slamch("Safe minimum");
231  eps = lapackf77_slamch("Precision");
232  smlnum = safmin / eps;
233  bignum = 1. / smlnum;
234  rmin = magma_ssqrt(smlnum);
235  rmax = magma_ssqrt(bignum);
236 
237  /* Scale matrix to allowable range, if necessary. */
238  anrm = lapackf77_slansy("M", uplo_, &n, a, &lda, work );
239  iscale = 0;
240  if (anrm > 0. && anrm < rmin) {
241  iscale = 1;
242  sigma = rmin / anrm;
243  } else if (anrm > rmax) {
244  iscale = 1;
245  sigma = rmax / anrm;
246  }
247  if (iscale == 1) {
248  lapackf77_slascl(uplo_, &izero, &izero, &d_one, &sigma, &n, &n, a,
249  &lda, info);
250  }
251 
252  /* Call SSYTRD to reduce symmetric matrix to tridiagonal form. */
253  // ssytrd work: e (n) + tau (n) + llwork (n*nb) ==> 2n + n*nb
254  // sstedx work: e (n) + tau (n) + z (n*n) + llwrk2 (1 + 4*n + n^2) ==> 1 + 6n + 2n^2
255  inde = 0;
256  indtau = inde + n;
257  indwrk = indtau + n;
258  indwk2 = indwrk + n*n;
259  llwork = lwork - indwrk;
260  llwrk2 = lwork - indwk2;
261 
262 //
263 #ifdef ENABLE_TIMER
264  magma_timestr_t start, end;
265  start = get_current_time();
266 #endif
267 
268  magma_ssytrd(uplo, n, a, lda, w, &work[inde],
269  &work[indtau], &work[indwrk], llwork, &iinfo);
270 
271 #ifdef ENABLE_TIMER
272  end = get_current_time();
273  printf("time ssytrd = %6.2f\n", GetTimerValue(start,end)/1000.);
274 #endif
275 
276  /* For eigenvalues only, call SSTERF. For eigenvectors, first call
277  SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
278  tridiagonal matrix, then call SORMTR to multiply it to the Householder
279  transformations represented as Householder vectors in A. */
280  if (! wantz) {
281  lapackf77_ssterf(&n, w, &work[inde], info);
282  } else {
283 
284 #ifdef ENABLE_TIMER
285  start = get_current_time();
286 #endif
287 
288  if (MAGMA_SUCCESS != magma_smalloc( &dwork, 3*n*(n/2 + 1) )) {
289  *info = MAGMA_ERR_DEVICE_ALLOC;
290  return *info;
291  }
292 
293  // TTT Possible bug for n < 128
294  magma_sstedx('A', n, 0., 0., 0, 0, w, &work[inde],
295  &work[indwrk], n, &work[indwk2],
296  llwrk2, iwork, liwork, dwork, info);
297 
298  magma_free( dwork );
299 
300 #ifdef ENABLE_TIMER
301  end = get_current_time();
302  printf("time sstedx = %6.2f\n", GetTimerValue(start,end)/1000.);
303  start = get_current_time();
304 #endif
305 
306  magma_sormtr(MagmaLeft, uplo, MagmaNoTrans, n, n, a, lda, &work[indtau],
307  &work[indwrk], n, &work[indwk2], llwrk2, &iinfo);
308 
309  lapackf77_slacpy("A", &n, &n, &work[indwrk], &n, a, &lda);
310 
311 #ifdef ENABLE_TIMER
312  end = get_current_time();
313  printf("time sormtr + copy = %6.2f\n", GetTimerValue(start,end)/1000.);
314 #endif
315  }
316 
317  /* If matrix was scaled, then rescale eigenvalues appropriately. */
318  if (iscale == 1) {
319  d__1 = 1. / sigma;
320  blasf77_sscal(&n, &d__1, w, &ione);
321  }
322 
323  work[0] = lwmin * (1. + lapackf77_slamch("Epsilon")); // round up
324  iwork[0] = liwmin;
325 
326  return *info;
327 } /* magma_ssyevd */
#define MagmaLeft
Definition: magma.h:68
#define __func__
Definition: common_magma.h:65
#define magma_ssqrt
Definition: common_magma.h:99
void blasf77_sscal(const magma_int_t *n, const float *alpha, float *x, const magma_int_t *incx)
#define MAGMA_ERR_DEVICE_ALLOC
Definition: magma_types.h:276
magma_int_t magma_sstedx(char range, magma_int_t n, float vl, float vu, magma_int_t il, magma_int_t iu, float *d, float *e, float *z, magma_int_t ldz, float *work, magma_int_t lwork, magma_int_t *iwork, magma_int_t liwork, float *dwork, magma_int_t *info)
Definition: sstedx.cpp:18
#define magma_free(ptr)
Definition: magma.h:57
void lapackf77_slacpy(const char *uplo, const magma_int_t *m, const magma_int_t *n, const float *A, const magma_int_t *lda, float *B, const magma_int_t *ldb)
int magma_int_t
Definition: magmablas.h:12
magma_int_t magma_ssytrd(char uplo, magma_int_t n, float *A, magma_int_t lda, float *d, float *e, float *tau, float *work, magma_int_t lwork, magma_int_t *info)
Definition: ssytrd.cpp:20
#define dwork(dev, i, j)
float lapackf77_slansy(const char *norm, const char *uplo, const magma_int_t *n, const float *A, const magma_int_t *lda, float *work)
void magma_xerbla(const char *srname, magma_int_t info)
Definition: xerbla.cpp:8
#define lapackf77_lsame
Definition: magma_lapack.h:23
#define MagmaLowerStr
Definition: magma.h:85
magma_int_t magma_get_ssytrd_nb(magma_int_t m)
Definition: get_nb.cpp:376
#define MAGMA_SUCCESS
Definition: magma.h:106
static magma_err_t magma_smalloc(magmaFloat_ptr *ptrPtr, size_t n)
Definition: magma.h:77
#define lapackf77_slamch
Definition: magma_lapack.h:26
#define MagmaNoVecStr
Definition: magma_types.h:423
double GetTimerValue(magma_timestr_t start, magma_timestr_t end)
Definition: timer.cpp:94
magma_int_t magma_sormtr(char side, char uplo, char trans, magma_int_t m, magma_int_t n, float *a, magma_int_t lda, float *tau, float *c, magma_int_t ldc, float *work, magma_int_t lwork, magma_int_t *info)
Definition: sormtr.cpp:17
#define MagmaUpperStr
Definition: magma.h:84
void lapackf77_ssyevd(const char *jobz, const char *uplo, const magma_int_t *n, float *A, const magma_int_t *lda, float *w, float *work, const magma_int_t *lwork, DWORKFORZ_AND_LD magma_int_t *iwork, const magma_int_t *liwork, magma_int_t *info)
#define MagmaNoTrans
Definition: magma.h:57
#define MagmaVecStr
Definition: magma_types.h:424
#define max(a, b)
Definition: common_magma.h:82
magma_timestr_t get_current_time(void)
Definition: timer.cpp:76
void lapackf77_slascl(const char *type, const magma_int_t *kl, const magma_int_t *ku, float *cfrom, float *cto, const magma_int_t *m, const magma_int_t *n, float *A, const magma_int_t *lda, magma_int_t *info)
#define lapackf77_ssterf
Definition: magma_lapack.h:34

Here is the call graph for this function:

Here is the caller graph for this function: