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

Go to the source code of this file.

Macros

#define PRECISION_c
 
#define VERSION3
 
#define vl(i, j)   (vl + (i) + (j)*ldvl)
 
#define vr(i, j)   (vr + (i) + (j)*ldvr)
 

Functions

magma_int_t magma_cgeev (char jobvl, char jobvr, magma_int_t n, magmaFloatComplex *A, magma_int_t lda, magmaFloatComplex *W, magmaFloatComplex *vl, magma_int_t ldvl, magmaFloatComplex *vr, magma_int_t ldvr, magmaFloatComplex *work, magma_int_t lwork, float *rwork, magma_int_t *info)
 

Macro Definition Documentation

#define PRECISION_c

Definition at line 15 of file cgeev.cpp.

#define VERSION3

Definition at line 22 of file cgeev.cpp.

#define vl (   i,
 
)    (vl + (i) + (j)*ldvl)
#define vr (   i,
 
)    (vr + (i) + (j)*ldvr)

Function Documentation

magma_int_t magma_cgeev ( char  jobvl,
char  jobvr,
magma_int_t  n,
magmaFloatComplex *  A,
magma_int_t  lda,
magmaFloatComplex *  W,
magmaFloatComplex *  vl,
magma_int_t  ldvl,
magmaFloatComplex *  vr,
magma_int_t  ldvr,
magmaFloatComplex *  work,
magma_int_t  lwork,
float *  rwork,
magma_int_t info 
)

Definition at line 25 of file cgeev.cpp.

References __func__, cblas_cscal(), cblas_csscal(), cblas_isamax(), CBLAS_SADDR, cblas_scnrm2(), dT, lapackf77_cgebak(), lapackf77_cgebal(), lapackf77_cgehrd(), lapackf77_chseqr(), lapackf77_clacpy(), lapackf77_clange(), lapackf77_clascl(), lapackf77_ctrevc(), lapackf77_cunghr(), lapackf77_lsame, lapackf77_slabad, lapackf77_slamch, MAGMA_C_CNJG, MAGMA_C_IMAG, MAGMA_C_MAKE, MAGMA_C_REAL, MAGMA_C_SET2REAL, magma_cgehrd(), magma_cgehrd2(), magma_cmalloc(), magma_cunghr(), MAGMA_ERR_DEVICE_ALLOC, magma_free, magma_get_cgehrd_nb(), magma_ssqrt, MAGMA_SUCCESS, magma_xerbla(), MagmaLowerStr, max, vl, and vr.

33 {
34 /* -- MAGMA (version 1.4.0) --
35  Univ. of Tennessee, Knoxville
36  Univ. of California, Berkeley
37  Univ. of Colorado, Denver
38  August 2013
39 
40  Purpose
41  =======
42  CGEEV computes for an N-by-N complex nonsymmetric matrix A, the
43  eigenvalues and, optionally, the left and/or right eigenvectors.
44 
45  The right eigenvector v(j) of A satisfies
46  A * v(j) = lambda(j) * v(j)
47  where lambda(j) is its eigenvalue.
48  The left eigenvector u(j) of A satisfies
49  u(j)**H * A = lambda(j) * u(j)**H
50  where u(j)**H denotes the conjugate transpose of u(j).
51 
52  The computed eigenvectors are normalized to have Euclidean norm
53  equal to 1 and largest component real.
54 
55  Arguments
56  =========
57  JOBVL (input) CHARACTER*1
58  = 'N': left eigenvectors of A are not computed;
59  = 'V': left eigenvectors of are computed.
60 
61  JOBVR (input) CHARACTER*1
62  = 'N': right eigenvectors of A are not computed;
63  = 'V': right eigenvectors of A are computed.
64 
65  N (input) INTEGER
66  The order of the matrix A. N >= 0.
67 
68  A (input/output) COMPLEX array, dimension (LDA,N)
69  On entry, the N-by-N matrix A.
70  On exit, A has been overwritten.
71 
72  LDA (input) INTEGER
73  The leading dimension of the array A. LDA >= max(1,N).
74 
75  W (output) COMPLEX array, dimension (N)
76  W contains the computed eigenvalues.
77 
78  VL (output) COMPLEX array, dimension (LDVL,N)
79  If JOBVL = 'V', the left eigenvectors u(j) are stored one
80  after another in the columns of VL, in the same order
81  as their eigenvalues.
82  If JOBVL = 'N', VL is not referenced.
83  u(j) = VL(:,j), the j-th column of VL.
84 
85  LDVL (input) INTEGER
86  The leading dimension of the array VL. LDVL >= 1; if
87  JOBVL = 'V', LDVL >= N.
88 
89  VR (output) COMPLEX array, dimension (LDVR,N)
90  If JOBVR = 'V', the right eigenvectors v(j) are stored one
91  after another in the columns of VR, in the same order
92  as their eigenvalues.
93  If JOBVR = 'N', VR is not referenced.
94  v(j) = VR(:,j), the j-th column of VR.
95 
96  LDVR (input) INTEGER
97  The leading dimension of the array VR. LDVR >= 1; if
98  JOBVR = 'V', LDVR >= N.
99 
100  WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
101  On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
102 
103  LWORK (input) INTEGER
104  The dimension of the array WORK. LWORK >= (1+nb)*N.
105 
106  If LWORK = -1, then a workspace query is assumed; the routine
107  only calculates the optimal size of the WORK array, returns
108  this value as the first entry of the WORK array, and no error
109  message related to LWORK is issued by XERBLA.
110 
111  RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
112 
113  INFO (output) INTEGER
114  = 0: successful exit
115  < 0: if INFO = -i, the i-th argument had an illegal value.
116  > 0: if INFO = i, the QR algorithm failed to compute all the
117  eigenvalues, and no eigenvectors have been computed;
118  elements and i+1:N of W contain eigenvalues which have
119  converged.
120  ===================================================================== */
121 
122  #define vl(i,j) (vl + (i) + (j)*ldvl)
123  #define vr(i,j) (vr + (i) + (j)*ldvr)
124 
125  magma_int_t c_one = 1;
126  magma_int_t c_zero = 0;
127 
128  float d__1, d__2;
129  magmaFloatComplex z__1, z__2;
130  magmaFloatComplex tmp;
131  float scl;
132  float dum[1], eps;
133  float anrm, cscale, bignum, smlnum;
134  magma_int_t i, k, ilo, ihi;
135  magma_int_t ibal, ierr, itau, iwrk, nout, liwrk, i__1, i__2, nb;
136  magma_int_t scalea, minwrk, irwork, lquery, wantvl, wantvr, select[1];
137 
138  char side[2] = {0, 0};
139  char jobvl_[2] = {jobvl, 0};
140  char jobvr_[2] = {jobvr, 0};
141 
142  irwork = 0;
143  *info = 0;
144  lquery = lwork == -1;
145  wantvl = lapackf77_lsame( jobvl_, "V" );
146  wantvr = lapackf77_lsame( jobvr_, "V" );
147  if (! wantvl && ! lapackf77_lsame( jobvl_, "N" )) {
148  *info = -1;
149  } else if (! wantvr && ! lapackf77_lsame( jobvr_, "N" )) {
150  *info = -2;
151  } else if (n < 0) {
152  *info = -3;
153  } else if (lda < max(1,n)) {
154  *info = -5;
155  } else if ( (ldvl < 1) || (wantvl && (ldvl < n))) {
156  *info = -8;
157  } else if ( (ldvr < 1) || (wantvr && (ldvr < n))) {
158  *info = -10;
159  }
160 
161  /* Compute workspace */
162  nb = magma_get_cgehrd_nb( n );
163  if (*info == 0) {
164  minwrk = (1+nb)*n;
165  work[0] = MAGMA_C_MAKE( (float) minwrk, 0. );
166 
167  if (lwork < minwrk && ! lquery) {
168  *info = -12;
169  }
170  }
171 
172  if (*info != 0) {
173  magma_xerbla( __func__, -(*info) );
174  return *info;
175  }
176  else if (lquery) {
177  return *info;
178  }
179 
180  /* Quick return if possible */
181  if (n == 0) {
182  return *info;
183  }
184 
185  #if defined(VERSION3)
186  magmaFloatComplex *dT;
187  if (MAGMA_SUCCESS != magma_cmalloc( &dT, nb*n )) {
188  *info = MAGMA_ERR_DEVICE_ALLOC;
189  return *info;
190  }
191  #endif
192 
193  /* Get machine constants */
194  eps = lapackf77_slamch( "P" );
195  smlnum = lapackf77_slamch( "S" );
196  bignum = 1. / smlnum;
197  lapackf77_slabad( &smlnum, &bignum );
198  smlnum = magma_ssqrt( smlnum ) / eps;
199  bignum = 1. / smlnum;
200 
201  /* Scale A if max element outside range [SMLNUM,BIGNUM] */
202  anrm = lapackf77_clange( "M", &n, &n, A, &lda, dum );
203  scalea = 0;
204  if (anrm > 0. && anrm < smlnum) {
205  scalea = 1;
206  cscale = smlnum;
207  } else if (anrm > bignum) {
208  scalea = 1;
209  cscale = bignum;
210  }
211  if (scalea) {
212  lapackf77_clascl( "G", &c_zero, &c_zero, &anrm, &cscale, &n, &n, A, &lda, &ierr );
213  }
214 
215  /* Balance the matrix
216  * (CWorkspace: none)
217  * (RWorkspace: need N) */
218  ibal = 0;
219  lapackf77_cgebal( "B", &n, A, &lda, &ilo, &ihi, &rwork[ibal], &ierr );
220 
221  /* Reduce to upper Hessenberg form
222  * (CWorkspace: need 2*N, prefer N + N*NB)
223  * (RWorkspace: none) */
224  itau = 0;
225  iwrk = itau + n;
226  liwrk = lwork - iwrk;
227 
228  #if defined(VERSION1)
229  // Version 1 - LAPACK
230  lapackf77_cgehrd( &n, &ilo, &ihi, A, &lda,
231  &work[itau], &work[iwrk], &liwrk, &ierr );
232  #elif defined(VERSION2)
233  // Version 2 - LAPACK consistent HRD
234  magma_cgehrd2( n, ilo, ihi, A, lda,
235  &work[itau], &work[iwrk], liwrk, &ierr );
236  #elif defined(VERSION3)
237  // Version 3 - LAPACK consistent MAGMA HRD + matrices T stored,
238  magma_cgehrd( n, ilo, ihi, A, lda,
239  &work[itau], &work[iwrk], liwrk, dT, &ierr );
240  #endif
241 
242  if (wantvl) {
243  /* Want left eigenvectors
244  * Copy Householder vectors to VL */
245  side[0] = 'L';
246  lapackf77_clacpy( MagmaLowerStr, &n, &n, A, &lda, vl, &ldvl );
247 
248  /* Generate unitary matrix in VL
249  * (CWorkspace: need 2*N-1, prefer N + (N-1)*NB)
250  * (RWorkspace: none) */
251  #if defined(VERSION1) || defined(VERSION2)
252  // Version 1 & 2 - LAPACK
253  lapackf77_cunghr( &n, &ilo, &ihi, vl, &ldvl, &work[itau],
254  &work[iwrk], &liwrk, &ierr );
255  #elif defined(VERSION3)
256  // Version 3 - LAPACK consistent MAGMA HRD + matrices T stored
257  magma_cunghr( n, ilo, ihi, vl, ldvl, &work[itau], dT, nb, &ierr );
258  #endif
259 
260  /* Perform QR iteration, accumulating Schur vectors in VL
261  * (CWorkspace: need 1, prefer HSWORK (see comments) )
262  * (RWorkspace: none) */
263  iwrk = itau;
264  liwrk = lwork - iwrk;
265  lapackf77_chseqr( "S", "V", &n, &ilo, &ihi, A, &lda, W,
266  vl, &ldvl, &work[iwrk], &liwrk, info );
267 
268  if (wantvr) {
269  /* Want left and right eigenvectors
270  * Copy Schur vectors to VR */
271  side[0] = 'B';
272  lapackf77_clacpy( "F", &n, &n, vl, &ldvl, vr, &ldvr );
273  }
274  }
275  else if (wantvr) {
276  /* Want right eigenvectors
277  * Copy Householder vectors to VR */
278  side[0] = 'R';
279  lapackf77_clacpy( "L", &n, &n, A, &lda, vr, &ldvr );
280 
281  /* Generate unitary matrix in VR
282  * (CWorkspace: need 2*N-1, prefer N + (N-1)*NB)
283  * (RWorkspace: none) */
284  #if defined(VERSION1) || defined(VERSION2)
285  // Version 1 & 2 - LAPACK
286  lapackf77_cunghr( &n, &ilo, &ihi, vr, &ldvr, &work[itau],
287  &work[iwrk], &liwrk, &ierr );
288  #elif defined(VERSION3)
289  // Version 3 - LAPACK consistent MAGMA HRD + matrices T stored
290  magma_cunghr( n, ilo, ihi, vr, ldvr, &work[itau], dT, nb, &ierr );
291  #endif
292 
293  /* Perform QR iteration, accumulating Schur vectors in VR
294  * (CWorkspace: need 1, prefer HSWORK (see comments) )
295  * (RWorkspace: none) */
296  iwrk = itau;
297  liwrk = lwork - iwrk;
298  lapackf77_chseqr( "S", "V", &n, &ilo, &ihi, A, &lda, W,
299  vr, &ldvr, &work[iwrk], &liwrk, info );
300  }
301  else {
302  /* Compute eigenvalues only
303  * (CWorkspace: need 1, prefer HSWORK (see comments) )
304  * (RWorkspace: none) */
305  iwrk = itau;
306  liwrk = lwork - iwrk;
307  lapackf77_chseqr( "E", "N", &n, &ilo, &ihi, A, &lda, W,
308  vr, &ldvr, &work[iwrk], &liwrk, info );
309  }
310 
311  /* If INFO > 0 from CHSEQR, then quit */
312  if (*info > 0) {
313  goto CLEANUP;
314  }
315 
316  if (wantvl || wantvr) {
317  /* Compute left and/or right eigenvectors
318  * (CWorkspace: need 2*N)
319  * (RWorkspace: need 2*N) */
320  irwork = ibal + n;
321  lapackf77_ctrevc( side, "B", select, &n, A, &lda, vl, &ldvl,
322  vr, &ldvr, &n, &nout, &work[iwrk], &rwork[irwork], &ierr );
323  }
324 
325  if (wantvl) {
326  /* Undo balancing of left eigenvectors
327  * (CWorkspace: none)
328  * (RWorkspace: need N) */
329  lapackf77_cgebak( "B", "L", &n, &ilo, &ihi, &rwork[ibal], &n,
330  vl, &ldvl, &ierr );
331 
332  /* Normalize left eigenvectors and make largest component real */
333  for (i = 0; i < n; ++i) {
334  scl = 1. / cblas_scnrm2( n, vl(0,i), 1 );
335  cblas_csscal( n, scl, vl(0,i), 1 );
336  for (k = 0; k < n; ++k) {
337  /* Computing 2nd power */
338  d__1 = MAGMA_C_REAL( *vl(k,i) );
339  d__2 = MAGMA_C_IMAG( *vl(k,i) );
340  rwork[irwork + k] = d__1*d__1 + d__2*d__2;
341  }
342  k = cblas_isamax( n, &rwork[irwork], 1 );
343  z__2 = MAGMA_C_CNJG( *vl(k,i) );
344  d__1 = magma_ssqrt( rwork[irwork + k] );
345  MAGMA_C_SSCALE( z__1, z__2, d__1 );
346  tmp = z__1;
347  cblas_cscal( n, CBLAS_SADDR(tmp), vl(0,i), 1 );
348  d__1 = MAGMA_C_REAL( *vl(k,i) );
349  MAGMA_C_SET2REAL( z__1, d__1 );
350  *vl(k,i) = z__1;
351  }
352  }
353 
354  if (wantvr) {
355  /* Undo balancing of right eigenvectors
356  * (CWorkspace: none)
357  * (RWorkspace: need N) */
358  lapackf77_cgebak( "B", "R", &n, &ilo, &ihi, &rwork[ibal], &n,
359  vr, &ldvr, &ierr );
360 
361  /* Normalize right eigenvectors and make largest component real */
362  for (i = 0; i < n; ++i) {
363  scl = 1. / cblas_scnrm2( n, vr(0,i), 1 );
364  cblas_csscal( n, scl, vr(0,i), 1 );
365  for (k = 0; k < n; ++k) {
366  /* Computing 2nd power */
367  d__1 = MAGMA_C_REAL( *vr(k,i) );
368  d__2 = MAGMA_C_IMAG( *vr(k,i) );
369  rwork[irwork + k] = d__1*d__1 + d__2*d__2;
370  }
371  k = cblas_isamax( n, &rwork[irwork], 1 );
372  z__2 = MAGMA_C_CNJG( *vr(k,i) );
373  d__1 = magma_ssqrt( rwork[irwork + k] );
374  MAGMA_C_SSCALE( z__1, z__2, d__1 );
375  tmp = z__1;
376  cblas_cscal( n, CBLAS_SADDR(tmp), vr(0,i), 1 );
377  d__1 = MAGMA_C_REAL( *vr(k,i) );
378  MAGMA_C_SET2REAL( z__1, d__1 );
379  *vr(k,i) = z__1;
380  }
381  }
382 
383 CLEANUP:
384  /* Undo scaling if necessary */
385  if (scalea) {
386  i__1 = n - (*info);
387  i__2 = max( n - (*info), 1 );
388  lapackf77_clascl( "G", &c_zero, &c_zero, &cscale, &anrm, &i__1, &c_one,
389  W + (*info), &i__2, &ierr );
390  if (*info > 0) {
391  i__1 = ilo - 1;
392  lapackf77_clascl( "G", &c_zero, &c_zero, &cscale, &anrm, &i__1, &c_one,
393  W, &n, &ierr );
394  }
395  }
396 
397  #if defined(VERSION3)
398  magma_free( dT );
399  #endif
400 
401  return *info;
402 } /* magma_cgeev */
float cblas_scnrm2(const int N, const void *X, const int incX)
void lapackf77_cgehrd(const magma_int_t *n, const magma_int_t *ilo, const magma_int_t *ihi, magmaFloatComplex *A, const magma_int_t *lda, magmaFloatComplex *tau, magmaFloatComplex *work, const magma_int_t *lwork, magma_int_t *info)
void lapackf77_ctrevc(const char *side, const char *howmny, magma_int_t *select, const magma_int_t *n, magmaFloatComplex *T, const magma_int_t *ldt, magmaFloatComplex *Vl, const magma_int_t *ldvl, magmaFloatComplex *Vr, const magma_int_t *ldvr, const magma_int_t *mm, magma_int_t *m, magmaFloatComplex *work, DWORKFORZ magma_int_t *info)
#define __func__
Definition: common_magma.h:65
magma_int_t magma_cgehrd2(magma_int_t n, magma_int_t ilo, magma_int_t ihi, magmaFloatComplex *A, magma_int_t lda, magmaFloatComplex *tau, magmaFloatComplex *work, magma_int_t lwork, magma_int_t *info)
Definition: cgehrd2.cpp:14
#define magma_ssqrt
Definition: common_magma.h:99
void lapackf77_chseqr(const char *job, const char *compz, const magma_int_t *n, const magma_int_t *ilo, const magma_int_t *ihi, magmaFloatComplex *H, const magma_int_t *ldh, WSPLIT, magmaFloatComplex *Z, const magma_int_t *ldz, magmaFloatComplex *work, const magma_int_t *lwork, magma_int_t *info)
static magma_err_t magma_cmalloc(magmaFloatComplex_ptr *ptrPtr, size_t n)
Definition: magma.h:79
#define MAGMA_ERR_DEVICE_ALLOC
Definition: magma_types.h:276
#define magma_free(ptr)
Definition: magma.h:57
#define MAGMA_C_IMAG(a)
Definition: magma.h:147
void cblas_csscal(const int N, const float alpha, void *X, const int incX)
int magma_int_t
Definition: magmablas.h:12
#define MAGMA_C_CNJG(v, t)
Definition: magma.h:142
#define W(k, n)
Definition: zgeqrf_mc.cpp:15
magma_int_t magma_get_cgehrd_nb(magma_int_t m)
Definition: get_nb.cpp:361
#define vl(i, j)
float lapackf77_clange(const char *norm, const magma_int_t *m, const magma_int_t *n, const magmaFloatComplex *A, const magma_int_t *lda, float *work)
void cblas_cscal(const int N, const void *alpha, void *X, const int incX)
magma_int_t magma_cunghr(magma_int_t n, magma_int_t ilo, magma_int_t ihi, magmaFloatComplex *a, magma_int_t lda, magmaFloatComplex *tau, magmaFloatComplex *dT, magma_int_t nb, magma_int_t *info)
Definition: cunghr.cpp:14
void magma_xerbla(const char *srname, magma_int_t info)
Definition: xerbla.cpp:8
#define lapackf77_lsame
Definition: magma_lapack.h:23
void lapackf77_cunghr(const magma_int_t *n, const magma_int_t *ilo, const magma_int_t *ihi, magmaFloatComplex *A, const magma_int_t *lda, const magmaFloatComplex *tau, magmaFloatComplex *work, const magma_int_t *lwork, magma_int_t *info)
void lapackf77_cgebal(const char *job, const magma_int_t *n, magmaFloatComplex *A, const magma_int_t *lda, magma_int_t *ilo, magma_int_t *ihi, float *scale, magma_int_t *info)
void lapackf77_clacpy(const char *uplo, const magma_int_t *m, const magma_int_t *n, const magmaFloatComplex *A, const magma_int_t *lda, magmaFloatComplex *B, const magma_int_t *ldb)
#define MagmaLowerStr
Definition: magma.h:85
#define lapackf77_slabad
Definition: magma_lapack.h:28
#define A(i, j)
Definition: cprint.cpp:16
void lapackf77_clascl(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, magmaFloatComplex *A, const magma_int_t *lda, magma_int_t *info)
#define MAGMA_SUCCESS
Definition: magma.h:106
#define lapackf77_slamch
Definition: magma_lapack.h:26
CBLAS_INDEX cblas_isamax(const int N, const float *X, const int incX)
#define MAGMA_C_REAL(a)
Definition: magma.h:146
magma_int_t magma_cgehrd(magma_int_t n, magma_int_t ilo, magma_int_t ihi, magmaFloatComplex *A, magma_int_t lda, magmaFloatComplex *tau, magmaFloatComplex *work, magma_int_t lwork, magmaFloatComplex *dT, magma_int_t *info)
Definition: cgehrd.cpp:17
void lapackf77_cgebak(const char *job, const char *side, const magma_int_t *n, const magma_int_t *ilo, const magma_int_t *ihi, const float *scale, const magma_int_t *m, magmaFloatComplex *V, const magma_int_t *ldv, magma_int_t *info)
#define dT(m)
#define vr(i, j)
#define MAGMA_C_SET2REAL(v, t)
Definition: magma.h:137
#define CBLAS_SADDR(a)
Definition: magma.h:204
#define max(a, b)
Definition: common_magma.h:82
#define MAGMA_C_MAKE(r, i)
Definition: magma.h:145

Here is the call graph for this function:

Here is the caller graph for this function: