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

Go to the source code of this file.

Macros

#define Z(ix, iy)   (z + (ix) + ldz * (iy))
 

Functions

magma_int_t magma_dstedx (char range, magma_int_t n, double vl, double vu, magma_int_t il, magma_int_t iu, double *d, double *e, double *z, magma_int_t ldz, double *work, magma_int_t lwork, magma_int_t *iwork, magma_int_t liwork, double *dwork, magma_int_t *info)
 

Macro Definition Documentation

#define Z (   ix,
  iy 
)    (z + (ix) + ldz * (iy))

Definition at line 14 of file dstedx.cpp.

Function Documentation

magma_int_t magma_dstedx ( char  range,
magma_int_t  n,
double  vl,
double  vu,
magma_int_t  il,
magma_int_t  iu,
double *  d,
double *  e,
double *  z,
magma_int_t  ldz,
double *  work,
magma_int_t  lwork,
magma_int_t iwork,
magma_int_t  liwork,
double *  dwork,
magma_int_t info 
)

Definition at line 18 of file dstedx.cpp.

References __func__, blasf77_dswap(), lapackf77_dlamch, lapackf77_dlanst(), lapackf77_dlascl(), lapackf77_dlaset(), lapackf77_dsteqr(), lapackf77_lsame, MAGMA_D_ABS, magma_dlaex0(), MAGMA_ERR_ILLEGAL_VALUE, magma_get_numthreads(), magma_get_smlsize_divideconquer(), magma_setlapack_numthreads(), MAGMA_SUCCESS, magma_xerbla(), max, min, and Z.

22 {
23 /* -- MAGMA (version 1.4.0) --
24  Univ. of Tennessee, Knoxville
25  Univ. of California, Berkeley
26  Univ. of Colorado, Denver
27  August 2013
28 
29  .. Scalar Arguments ..
30  CHARACTER RANGE
31  INTEGER IL, IU, INFO, LDZ, LIWORK, LWORK, N
32  DOUBLE PRECISION VL, VU
33  ..
34  .. Array Arguments ..
35  INTEGER IWORK( * )
36  DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * ),
37  $ DWORK ( * )
38  ..
39 
40  Purpose
41  =======
42  DSTEDX computes some eigenvalues and, optionally, eigenvectors of a
43  symmetric tridiagonal matrix using the divide and conquer method.
44 
45  This code makes very mild assumptions about floating point
46  arithmetic. It will work on machines with a guard digit in
47  add/subtract, or on those binary machines without guard digits
48  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
49  It could conceivably fail on hexadecimal or decimal machines
50  without guard digits, but we know of none. See DLAEX3 for details.
51 
52  Arguments
53  =========
54  RANGE (input) CHARACTER*1
55  = 'A': all eigenvalues will be found.
56  = 'V': all eigenvalues in the half-open interval (VL,VU]
57  will be found.
58  = 'I': the IL-th through IU-th eigenvalues will be found.
59 
60  N (input) INTEGER
61  The dimension of the symmetric tridiagonal matrix. N >= 0.
62 
63  VL (input) DOUBLE PRECISION
64  VU (input) DOUBLE PRECISION
65  If RANGE='V', the lower and upper bounds of the interval to
66  be searched for eigenvalues. VL < VU.
67  Not referenced if RANGE = 'A' or 'I'.
68 
69  IL (input) INTEGER
70  IU (input) INTEGER
71  If RANGE='I', the indices (in ascending order) of the
72  smallest and largest eigenvalues to be returned.
73  1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
74  Not referenced if RANGE = 'A' or 'V'.
75 
76  D (input/output) DOUBLE PRECISION array, dimension (N)
77  On entry, the diagonal elements of the tridiagonal matrix.
78  On exit, if INFO = 0, the eigenvalues in ascending order.
79 
80  E (input/output) DOUBLE PRECISION array, dimension (N-1)
81  On entry, the subdiagonal elements of the tridiagonal matrix.
82  On exit, E has been destroyed.
83 
84  Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
85  On exit, if INFO = 0, Z contains the orthonormal eigenvectors
86  of the symmetric tridiagonal matrix.
87 
88  LDZ (input) INTEGER
89  The leading dimension of the array Z. LDZ >= max(1,N).
90 
91  WORK (workspace/output) DOUBLE PRECISION array,
92  dimension (LWORK)
93  On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
94 
95  LWORK (input) INTEGER
96  The dimension of the array WORK.
97  If N > 1 then LWORK >= ( 1 + 4*N + N**2 ).
98  Note that if N is less than or
99  equal to the minimum divide size, usually 25, then LWORK need
100  only be max(1,2*(N-1)).
101 
102  If LWORK = -1, then a workspace query is assumed; the routine
103  only calculates the optimal size of the WORK array, returns
104  this value as the first entry of the WORK array, and no error
105  message related to LWORK is issued by XERBLA.
106 
107  IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
108  On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
109 
110  LIWORK (input) INTEGER
111  The dimension of the array IWORK.
112  LIWORK >= ( 3 + 5*N ).
113  Note that if N is less than or
114  equal to the minimum divide size, usually 25, then LIWORK
115  need only be 1.
116 
117  If LIWORK = -1, then a workspace query is assumed; the
118  routine only calculates the optimal size of the IWORK array,
119  returns this value as the first entry of the IWORK array, and
120  no error message related to LIWORK is issued by XERBLA.
121 
122  DWORK (device workspace) DOUBLE PRECISION array, dimension (3*N*N/2+3*N)
123 
124  INFO (output) INTEGER
125  = 0: successful exit.
126  < 0: if INFO = -i, the i-th argument had an illegal value.
127  > 0: The algorithm failed to compute an eigenvalue while
128  working on the submatrix lying in rows and columns
129  INFO/(N+1) through mod(INFO,N+1).
130 
131  Further Details
132  ===============
133  Based on contributions by
134  Jeff Rutter, Computer Science Division, University of California
135  at Berkeley, USA
136  Modified by Francoise Tisseur, University of Tennessee.
137 
138  =====================================================================
139 */
140  char range_[2] = {range, 0};
141 
142  double d_zero = 0.;
143  double d_one = 1.;
144  magma_int_t izero = 0;
145  magma_int_t ione = 1;
146 
147 
148  magma_int_t alleig, indeig, valeig, lquery;
149  magma_int_t i, j, k, m;
150  magma_int_t liwmin, lwmin;
151  magma_int_t start, end, smlsiz;
152  double eps, orgnrm, p, tiny;
153 
154  // Test the input parameters.
155 
156  alleig = lapackf77_lsame(range_, "A");
157  valeig = lapackf77_lsame(range_, "V");
158  indeig = lapackf77_lsame(range_, "I");
159  lquery = lwork == -1 || liwork == -1;
160 
161  *info = 0;
162 
163  if (! (alleig || valeig || indeig)) {
164  *info = -1;
165  } else if (n < 0) {
166  *info = -2;
167  } else if (ldz < max(1,n)) {
168  *info = -10;
169  } else {
170  if (valeig) {
171  if (n > 0 && vu <= vl) {
172  *info = -4;
173  }
174  } else if (indeig) {
175  if (il < 1 || il > max(1,n)) {
176  *info = -5;
177  } else if (iu < min(n,il) || iu > n) {
178  *info = -6;
179  }
180  }
181  }
182 
183  if (*info == 0) {
184  // Compute the workspace requirements
185 
187  if( n <= 1 ){
188  lwmin = 1;
189  liwmin = 1;
190  } else {
191  lwmin = 1 + 4*n + n*n;
192  liwmin = 3 + 5*n;
193  }
194 
195  work[0] = lwmin;
196  iwork[0] = liwmin;
197 
198  if (lwork < lwmin && ! lquery) {
199  *info = -12;
200  } else if (liwork < liwmin && ! lquery) {
201  *info = -14;
202  }
203  }
204 
205  if (*info != 0) {
206  magma_xerbla( __func__, -(*info));
208  } else if (lquery) {
209  return MAGMA_SUCCESS;
210  }
211 
212  // Quick return if possible
213 
214  if(n==0)
215  return MAGMA_SUCCESS;
216  if(n==1){
217  *z = 1.;
218  return MAGMA_SUCCESS;
219  }
220 
221  /* determine the number of threads */
222  magma_int_t threads = magma_get_numthreads();
224 
225 #ifdef ENABLE_DEBUG
226  printf(" D&C is using %d threads\n", threads);
227 #endif
228 
229  // If N is smaller than the minimum divide size (SMLSIZ+1), then
230  // solve the problem with another solver.
231 
232  if (n < smlsiz){
233  char char_I[]= {'I', 0};
234  lapackf77_dsteqr(char_I, &n, d, e, z, &ldz, work, info);
235  } else {
236  char char_F[]= {'F', 0};
237  lapackf77_dlaset(char_F, &n, &n, &d_zero, &d_one, z, &ldz);
238 
239  //Scale.
240  char char_M[]= {'M', 0};
241 
242  orgnrm = lapackf77_dlanst(char_M, &n, d, e);
243 
244  if (orgnrm == 0){
245  work[0] = lwmin;
246  iwork[0] = liwmin;
247  return MAGMA_SUCCESS;
248  }
249 
250  eps = lapackf77_dlamch( "Epsilon" );
251 
252  if (alleig){
253  start = 0;
254  while ( start < n ){
255 
256  // Let FINISH be the position of the next subdiagonal entry
257  // such that E( END ) <= TINY or FINISH = N if no such
258  // subdiagonal exists. The matrix identified by the elements
259  // between START and END constitutes an independent
260  // sub-problem.
261 
262  for(end = start+1; end < n; ++end){
263  tiny = eps * sqrt( MAGMA_D_ABS(d[end-1]*d[end]));
264  if (MAGMA_D_ABS(e[end-1]) <= tiny)
265  break;
266  }
267 
268  // (Sub) Problem determined. Compute its size and solve it.
269 
270  m = end - start;
271  if (m==1){
272  start = end;
273  continue;
274  }
275  if (m > smlsiz){
276 
277  // Scale
278  char char_G[] = {'G', 0};
279  orgnrm = lapackf77_dlanst(char_M, &m, &d[start], &e[start]);
280  lapackf77_dlascl(char_G, &izero, &izero, &orgnrm, &d_one, &m, &ione, &d[start], &m, info);
281  magma_int_t mm = m-1;
282  lapackf77_dlascl(char_G, &izero, &izero, &orgnrm, &d_one, &mm, &ione, &e[start], &mm, info);
283 
284  magma_dlaex0( m, &d[start], &e[start], Z(start, start), ldz, work, iwork, dwork, 'A', vl, vu, il, iu, info);
285 
286  if( *info != 0) {
287  return MAGMA_SUCCESS;
288  }
289 
290  // Scale Back
291  lapackf77_dlascl(char_G, &izero, &izero, &d_one, &orgnrm, &m, &ione, &d[start], &m, info);
292 
293  } else {
294 
295  char char_I[]= {'I', 0};
296  lapackf77_dsteqr( char_I, &m, &d[start], &e[start], Z(start, start), &ldz, work, info);
297  if (*info != 0){
298  *info = (start+1) *(n+1) + end;
299  }
300  }
301 
302  start = end;
303  }
304 
305 
306  // If the problem split any number of times, then the eigenvalues
307  // will not be properly ordered. Here we permute the eigenvalues
308  // (and the associated eigenvectors) into ascending order.
309 
310  if (m < n){
311 
312  // Use Selection Sort to minimize swaps of eigenvectors
313  for (i = 1; i < n; ++i){
314  k = i-1;
315  p = d[i-1];
316  for (j = i; j < n; ++j){
317  if (d[j] < p){
318  k = j;
319  p = d[j];
320  }
321  }
322  if(k != i-1) {
323  d[k] = d[i-1];
324  d[i-1] = p;
325  blasf77_dswap(&n, Z(0,i-1), &ione, Z(0,k), &ione);
326  }
327  }
328  }
329 
330  } else {
331 
332  // Scale
333  char char_G[] = {'G', 0};
334  lapackf77_dlascl(char_G, &izero, &izero, &orgnrm, &d_one, &n, &ione, d, &n, info);
335  magma_int_t nm = n-1;
336  lapackf77_dlascl(char_G, &izero, &izero, &orgnrm, &d_one, &nm, &ione, e, &nm, info);
337 
338  magma_dlaex0( n, d, e, z, ldz, work, iwork, dwork, range, vl, vu, il, iu, info);
339 
340  if( *info != 0) {
341  return MAGMA_SUCCESS;
342  }
343 
344  // Scale Back
345  lapackf77_dlascl(char_G, &izero, &izero, &d_one, &orgnrm, &n, &ione, d, &n, info);
346 
347  }
348  }
349 
350  work[0] = lwmin;
351  iwork[0] = liwmin;
352 
353  return MAGMA_SUCCESS;
354 
355 } /* dstedx */
#define MAGMA_ERR_ILLEGAL_VALUE
Definition: magma.h:107
#define min(a, b)
Definition: common_magma.h:86
#define __func__
Definition: common_magma.h:65
int magma_int_t
Definition: magmablas.h:12
#define lapackf77_dlamch
Definition: magma_lapack.h:27
#define dwork(dev, i, j)
#define vl(i, j)
double lapackf77_dlanst(const char *norm, const magma_int_t *n, const double *d, const double *e)
#define Z(ix, iy)
Definition: dstedx.cpp:14
void magma_setlapack_numthreads(magma_int_t num_threads)
void magma_xerbla(const char *srname, magma_int_t info)
Definition: xerbla.cpp:8
#define lapackf77_lsame
Definition: magma_lapack.h:23
void lapackf77_dlascl(const char *type, const magma_int_t *kl, const magma_int_t *ku, double *cfrom, double *cto, const magma_int_t *m, const magma_int_t *n, double *A, const magma_int_t *lda, magma_int_t *info)
#define MAGMA_SUCCESS
Definition: magma.h:106
void blasf77_dswap(const magma_int_t *n, double *x, const magma_int_t *incx, double *y, const magma_int_t *incy)
magma_int_t magma_dlaex0(magma_int_t n, double *d, double *e, double *q, magma_int_t ldq, double *work, magma_int_t *iwork, double *dwork, char range, double vl, double vu, magma_int_t il, magma_int_t iu, magma_int_t *info)
Definition: dlaex0.cpp:17
void lapackf77_dsteqr(const char *compz, const magma_int_t *n, double *d, double *e, double *Z, const magma_int_t *ldz, double *work, magma_int_t *info)
magma_int_t magma_get_smlsize_divideconquer()
Definition: get_nb.cpp:768
magma_int_t magma_get_numthreads()
void lapackf77_dlaset(const char *uplo, const magma_int_t *m, const magma_int_t *n, const double *alpha, const double *beta, double *A, const magma_int_t *lda)
#define max(a, b)
Definition: common_magma.h:82
#define MAGMA_D_ABS(a)
Definition: magma.h:174

Here is the call graph for this function:

Here is the caller graph for this function: