MAGMA  magma-1.4.0
Matrix Algebra on GPU and Multicore Architectures
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
dlaex0.cpp
Go to the documentation of this file.
1 /*
2  -- MAGMA (version 1.4.0) --
3  Univ. of Tennessee, Knoxville
4  Univ. of California, Berkeley
5  Univ. of Colorado, Denver
6  August 2013
7 
8  @author Raffaele Solca
9 
10  @precisions normal d -> s
11 */
12 #include "common_magma.h"
13 
14 #define Q(ix, iy) (q + (ix) + ldq * (iy))
15 
16 extern "C" magma_int_t
17 magma_dlaex0(magma_int_t n, double* d, double* e, double* q, magma_int_t ldq,
18  double* work, magma_int_t* iwork, double* dwork,
19  char range, double vl, double vu,
20  magma_int_t il, magma_int_t iu, magma_int_t* info)
21 {
22 /* -- MAGMA (version 1.4.0) --
23  Univ. of Tennessee, Knoxville
24  Univ. of California, Berkeley
25  Univ. of Colorado, Denver
26  August 2013
27 
28  .. Scalar Arguments ..
29  CHARACTER RANGE
30  INTEGER IL, IU, INFO, LDQ, N
31  DOUBLE PRECISION VL, VU
32  ..
33  .. Array Arguments ..
34  INTEGER IWORK( * )
35  DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ),
36  $ WORK( * ), DWORK( * )
37  ..
38 
39  Purpose
40  =======
41  DLAEX0 computes all eigenvalues and the choosen eigenvectors of a
42  symmetric tridiagonal matrix using the divide and conquer method.
43 
44  Arguments
45  =========
46  N (input) INTEGER
47  The dimension of the symmetric tridiagonal matrix. N >= 0.
48 
49  D (input/output) DOUBLE PRECISION array, dimension (N)
50  On entry, the main diagonal of the tridiagonal matrix.
51  On exit, its eigenvalues.
52 
53  E (input) DOUBLE PRECISION array, dimension (N-1)
54  The off-diagonal elements of the tridiagonal matrix.
55  On exit, E has been destroyed.
56 
57  Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
58  On entry, Q will be the identity matrix.
59  On exit, Q contains the eigenvectors of the
60  tridiagonal matrix.
61 
62  LDQ (input) INTEGER
63  The leading dimension of the array Q. If eigenvectors are
64  desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
65 
66  WORK (workspace) DOUBLE PRECISION array,
67  the dimension of WORK >= 4*N + N**2.
68 
69  IWORK (workspace) INTEGER array,
70  the dimension of IWORK >= 3 + 5*N.
71 
72  DWORK (device workspace) DOUBLE PRECISION array, dimension (3*N*N/2+3*N)
73 
74  RANGE (input) CHARACTER*1
75  = 'A': all eigenvalues will be found.
76  = 'V': all eigenvalues in the half-open interval (VL,VU]
77  will be found.
78  = 'I': the IL-th through IU-th eigenvalues will be found.
79 
80  VL (input) DOUBLE PRECISION
81  VU (input) DOUBLE PRECISION
82  If RANGE='V', the lower and upper bounds of the interval to
83  be searched for eigenvalues. VL < VU.
84  Not referenced if RANGE = 'A' or 'I'.
85 
86  IL (input) INTEGER
87  IU (input) INTEGER
88  If RANGE='I', the indices (in ascending order) of the
89  smallest and largest eigenvalues to be returned.
90  1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
91  Not referenced if RANGE = 'A' or 'V'.
92 
93  INFO (output) INTEGER
94  = 0: successful exit.
95  < 0: if INFO = -i, the i-th argument had an illegal value.
96  > 0: The algorithm failed to compute an eigenvalue while
97  working on the submatrix lying in rows and columns
98  INFO/(N+1) through mod(INFO,N+1).
99 
100  Further Details
101  ===============
102  Based on contributions by
103  Jeff Rutter, Computer Science Division, University of California
104  at Berkeley, USA
105 
106  ===================================================================== */
107 
108  magma_int_t ione = 1;
109  char range_;
110  magma_int_t curlvl, i, indxq;
111  magma_int_t j, k, matsiz, msd2, smlsiz;
112  magma_int_t submat, subpbs, tlvls;
113 
114 
115  // Test the input parameters.
116 
117  *info = 0;
118 
119  if( n < 0 )
120  *info = -1;
121  else if( ldq < max(1, n) )
122  *info = -5;
123  if( *info != 0 ){
124  magma_xerbla( __func__, -*info );
126  }
127 
128  // Quick return if possible
129  if(n == 0)
130  return MAGMA_SUCCESS;
131 
133 
134  // Determine the size and placement of the submatrices, and save in
135  // the leading elements of IWORK.
136 
137  iwork[0] = n;
138  subpbs= 1;
139  tlvls = 0;
140  while (iwork[subpbs - 1] > smlsiz) {
141  for (j = subpbs; j > 0; --j){
142  iwork[2*j - 1] = (iwork[j-1]+1)/2;
143  iwork[2*j - 2] = iwork[j-1]/2;
144  }
145  ++tlvls;
146  subpbs *= 2;
147  }
148  for (j=1; j<subpbs; ++j)
149  iwork[j] += iwork[j-1];
150 
151  // Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
152  // using rank-1 modifications (cuts).
153 
154  for(i=0; i < subpbs-1; ++i){
155  submat = iwork[i];
156  d[submat-1] -= MAGMA_D_ABS(e[submat-1]);
157  d[submat] -= MAGMA_D_ABS(e[submat-1]);
158  }
159 
160  indxq = 4*n + 3;
161 
162  // Solve each submatrix eigenproblem at the bottom of the divide and
163  // conquer tree.
164 
165  char char_I[] = {'I', 0};
166 
167 #ifdef ENABLE_TIMER_DIVIDE_AND_CONQUER
168  magma_timestr_t start, end;
169  start = get_current_time();
170 #endif
171 
172  for (i = 0; i < subpbs; ++i){
173  if(i == 0){
174  submat = 0;
175  matsiz = iwork[0];
176  } else {
177  submat = iwork[i-1];
178  matsiz = iwork[i] - iwork[i-1];
179  }
180  lapackf77_dsteqr(char_I , &matsiz, &d[submat], &e[submat],
181  Q(submat, submat), &ldq, work, info); // change to edc?
182  if(*info != 0){
183  printf("info: %d\n, submat: %d\n", (int) *info, (int) submat);
184  *info = (submat+1)*(n+1) + submat + matsiz;
185  printf("info: %d\n", (int) *info);
186  return MAGMA_SUCCESS;
187  }
188  k = 1;
189  for(j = submat; j < iwork[i]; ++j){
190  iwork[indxq+j] = k;
191  ++k;
192  }
193  }
194 
195 #ifdef ENABLE_TIMER_DIVIDE_AND_CONQUER
196  end = get_current_time();
197  printf(" for: dsteqr = %6.2f\n", GetTimerValue(start,end)/1000.);
198 #endif
199  // Successively merge eigensystems of adjacent submatrices
200  // into eigensystem for the corresponding larger matrix.
201 
202  curlvl = 1;
203  while (subpbs > 1){
204 #ifdef ENABLE_TIMER_DIVIDE_AND_CONQUER
205  magma_timestr_t start, end;
206  start = get_current_time();
207 #endif
208  for (i=0; i<subpbs-1; i+=2){
209  if(i == 0){
210  submat = 0;
211  matsiz = iwork[1];
212  msd2 = iwork[0];
213  } else {
214  submat = iwork[i-1];
215  matsiz = iwork[i+1] - iwork[i-1];
216  msd2 = matsiz / 2;
217  }
218 
219  // Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
220  // into an eigensystem of size MATSIZ.
221  // DLAEX1 is used only for the full eigensystem of a tridiagonal
222  // matrix.
223 
224  if (matsiz == n)
225  range_=range;
226  else
227  // We need all the eigenvectors if it is not last step
228  range_='A';
229 
230  magma_dlaex1(matsiz, &d[submat], Q(submat, submat), ldq,
231  &iwork[indxq+submat], e[submat+msd2-1], msd2,
232  work, &iwork[subpbs], dwork,
233  range_, vl, vu, il, iu, info);
234 
235  if(*info != 0){
236  *info = (submat+1)*(n+1) + submat + matsiz;
237  return MAGMA_SUCCESS;
238  }
239  iwork[i/2]= iwork[i+1];
240  }
241  subpbs /= 2;
242  ++curlvl;
243 #ifdef ENABLE_TIMER_DIVIDE_AND_CONQUER
244  end = get_current_time();
245  //printf("%d: time: %6.2f\n", curlvl, GetTimerValue(start,end)/1000.);
246 #endif
247 
248  }
249 
250  // Re-merge the eigenvalues/vectors which were deflated at the final
251  // merge step.
252 
253  for(i = 0; i<n; ++i){
254  j = iwork[indxq+i] - 1;
255  work[i] = d[j];
256  blasf77_dcopy(&n, Q(0, j), &ione, &work[ n*(i+1) ], &ione);
257  }
258  blasf77_dcopy(&n, work, &ione, d, &ione);
259  char char_A[] = {'A',0};
260  lapackf77_dlacpy ( char_A, &n, &n, &work[n], &n, q, &ldq );
261 
262  return MAGMA_SUCCESS;
263 
264 } /* magma_dlaex0 */
265 
#define MAGMA_ERR_ILLEGAL_VALUE
Definition: magma.h:107
void lapackf77_dlacpy(const char *uplo, const magma_int_t *m, const magma_int_t *n, const double *A, const magma_int_t *lda, double *B, const magma_int_t *ldb)
#define __func__
Definition: common_magma.h:65
int magma_int_t
Definition: magmablas.h:12
void blasf77_dcopy(const magma_int_t *n, const double *x, const magma_int_t *incx, double *y, const magma_int_t *incy)
#define dwork(dev, i, j)
#define vl(i, j)
magma_int_t magma_dlaex1(magma_int_t n, double *d, double *q, magma_int_t ldq, magma_int_t *indxq, double rho, magma_int_t cutpnt, 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: dlaex1.cpp:17
void magma_xerbla(const char *srname, magma_int_t info)
Definition: xerbla.cpp:8
#define MAGMA_SUCCESS
Definition: magma.h:106
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
double GetTimerValue(magma_timestr_t start, magma_timestr_t end)
Definition: timer.cpp:94
#define max(a, b)
Definition: common_magma.h:82
#define MAGMA_D_ABS(a)
Definition: magma.h:174
magma_timestr_t get_current_time(void)
Definition: timer.cpp:76
#define Q(ix, iy)
Definition: dlaex0.cpp:14