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

Go to the source code of this file.

Macros

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

Functions

magma_int_t magma_zgeev_m (char jobvl, char jobvr, magma_int_t n, magmaDoubleComplex *A, magma_int_t lda, magmaDoubleComplex *W, magmaDoubleComplex *vl, magma_int_t ldvl, magmaDoubleComplex *vr, magma_int_t ldvr, magmaDoubleComplex *work, magma_int_t lwork, double *rwork, magma_int_t *info)
 

Macro Definition Documentation

#define PRECISION_z

Definition at line 15 of file zgeev_m.cpp.

#define Version5

Definition at line 24 of file zgeev_m.cpp.

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

Function Documentation

magma_int_t magma_zgeev_m ( char  jobvl,
char  jobvr,
magma_int_t  n,
magmaDoubleComplex *  A,
magma_int_t  lda,
magmaDoubleComplex *  W,
magmaDoubleComplex *  vl,
magma_int_t  ldvl,
magmaDoubleComplex *  vr,
magma_int_t  ldvr,
magmaDoubleComplex *  work,
magma_int_t  lwork,
double *  rwork,
magma_int_t info 
)

Definition at line 27 of file zgeev_m.cpp.

References __func__, cblas_dznrm2(), cblas_idamax(), CBLAS_SADDR, cblas_zdscal(), cblas_zscal(), dT, lapackf77_dlabad, lapackf77_dlamch, lapackf77_lsame, lapackf77_zgebak, lapackf77_zgebal, lapackf77_zgehrd, lapackf77_zhseqr, lapackf77_zlacpy, lapackf77_zlange, lapackf77_zlascl, lapackf77_ztrevc, lapackf77_zunghr, magma_dsqrt, MAGMA_ERR_DEVICE_ALLOC, MAGMA_ERR_HOST_ALLOC, magma_free, magma_free_cpu(), magma_get_zgehrd_nb(), MAGMA_SUCCESS, magma_xerbla(), MAGMA_Z_CNJG, MAGMA_Z_DSCALE, MAGMA_Z_IMAG, MAGMA_Z_MAKE, MAGMA_Z_REAL, MAGMA_Z_SET2REAL, magma_zgehrd(), magma_zgehrd2(), magma_zgehrd_m(), magma_zmalloc(), magma_zmalloc_cpu(), magma_zsetmatrix, magma_zunghr(), magma_zunghr_m(), MagmaLowerStr, max, T, vl, and vr.

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

Here is the call graph for this function:

Here is the caller graph for this function: