ergo
template_lapack_laebz.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_LAEBZ_HEADER
36 #define TEMPLATE_LAPACK_LAEBZ_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_laebz(const integer *ijob, const integer *nitmax, const integer *n,
41  const integer *mmax, const integer *minp, const integer *nbmin, const Treal *abstol,
42  const Treal *reltol, const Treal *pivmin, const Treal *d__, const Treal *
43  e, const Treal *e2, integer *nval, Treal *ab, Treal *c__,
44  integer *mout, integer *nab, Treal *work, integer *iwork,
45  integer *info)
46 {
47 /* -- LAPACK auxiliary routine (version 3.0) --
48  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
49  Courant Institute, Argonne National Lab, and Rice University
50  June 30, 1999
51 
52 
53  Purpose
54  =======
55 
56  DLAEBZ contains the iteration loops which compute and use the
57  function N(w), which is the count of eigenvalues of a symmetric
58  tridiagonal matrix T less than or equal to its argument w. It
59  performs a choice of two types of loops:
60 
61  IJOB=1, followed by
62  IJOB=2: It takes as input a list of intervals and returns a list of
63  sufficiently small intervals whose union contains the same
64  eigenvalues as the union of the original intervals.
65  The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
66  The output interval (AB(j,1),AB(j,2)] will contain
67  eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
68 
69  IJOB=3: It performs a binary search in each input interval
70  (AB(j,1),AB(j,2)] for a point w(j) such that
71  N(w(j))=NVAL(j), and uses C(j) as the starting point of
72  the search. If such a w(j) is found, then on output
73  AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output
74  (AB(j,1),AB(j,2)] will be a small interval containing the
75  point where N(w) jumps through NVAL(j), unless that point
76  lies outside the initial interval.
77 
78  Note that the intervals are in all cases half-open intervals,
79  i.e., of the form (a,b] , which includes b but not a .
80 
81  To avoid underflow, the matrix should be scaled so that its largest
82  element is no greater than overflow**(1/2) * underflow**(1/4)
83  in absolute value. To assure the most accurate computation
84  of small eigenvalues, the matrix should be scaled to be
85  not much smaller than that, either.
86 
87  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
88  Matrix", Report CS41, Computer Science Dept., Stanford
89  University, July 21, 1966
90 
91  Note: the arguments are, in general, *not* checked for unreasonable
92  values.
93 
94  Arguments
95  =========
96 
97  IJOB (input) INTEGER
98  Specifies what is to be done:
99  = 1: Compute NAB for the initial intervals.
100  = 2: Perform bisection iteration to find eigenvalues of T.
101  = 3: Perform bisection iteration to invert N(w), i.e.,
102  to find a point which has a specified number of
103  eigenvalues of T to its left.
104  Other values will cause DLAEBZ to return with INFO=-1.
105 
106  NITMAX (input) INTEGER
107  The maximum number of "levels" of bisection to be
108  performed, i.e., an interval of width W will not be made
109  smaller than 2^(-NITMAX) * W. If not all intervals
110  have converged after NITMAX iterations, then INFO is set
111  to the number of non-converged intervals.
112 
113  N (input) INTEGER
114  The dimension n of the tridiagonal matrix T. It must be at
115  least 1.
116 
117  MMAX (input) INTEGER
118  The maximum number of intervals. If more than MMAX intervals
119  are generated, then DLAEBZ will quit with INFO=MMAX+1.
120 
121  MINP (input) INTEGER
122  The initial number of intervals. It may not be greater than
123  MMAX.
124 
125  NBMIN (input) INTEGER
126  The smallest number of intervals that should be processed
127  using a vector loop. If zero, then only the scalar loop
128  will be used.
129 
130  ABSTOL (input) DOUBLE PRECISION
131  The minimum (absolute) width of an interval. When an
132  interval is narrower than ABSTOL, or than RELTOL times the
133  larger (in magnitude) endpoint, then it is considered to be
134  sufficiently small, i.e., converged. This must be at least
135  zero.
136 
137  RELTOL (input) DOUBLE PRECISION
138  The minimum relative width of an interval. When an interval
139  is narrower than ABSTOL, or than RELTOL times the larger (in
140  magnitude) endpoint, then it is considered to be
141  sufficiently small, i.e., converged. Note: this should
142  always be at least radix*machine epsilon.
143 
144  PIVMIN (input) DOUBLE PRECISION
145  The minimum absolute value of a "pivot" in the Sturm
146  sequence loop. This *must* be at least max |e(j)**2| *
147  safe_min and at least safe_min, where safe_min is at least
148  the smallest number that can divide one without overflow.
149 
150  D (input) DOUBLE PRECISION array, dimension (N)
151  The diagonal elements of the tridiagonal matrix T.
152 
153  E (input) DOUBLE PRECISION array, dimension (N)
154  The offdiagonal elements of the tridiagonal matrix T in
155  positions 1 through N-1. E(N) is arbitrary.
156 
157  E2 (input) DOUBLE PRECISION array, dimension (N)
158  The squares of the offdiagonal elements of the tridiagonal
159  matrix T. E2(N) is ignored.
160 
161  NVAL (input/output) INTEGER array, dimension (MINP)
162  If IJOB=1 or 2, not referenced.
163  If IJOB=3, the desired values of N(w). The elements of NVAL
164  will be reordered to correspond with the intervals in AB.
165  Thus, NVAL(j) on output will not, in general be the same as
166  NVAL(j) on input, but it will correspond with the interval
167  (AB(j,1),AB(j,2)] on output.
168 
169  AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2)
170  The endpoints of the intervals. AB(j,1) is a(j), the left
171  endpoint of the j-th interval, and AB(j,2) is b(j), the
172  right endpoint of the j-th interval. The input intervals
173  will, in general, be modified, split, and reordered by the
174  calculation.
175 
176  C (input/output) DOUBLE PRECISION array, dimension (MMAX)
177  If IJOB=1, ignored.
178  If IJOB=2, workspace.
179  If IJOB=3, then on input C(j) should be initialized to the
180  first search point in the binary search.
181 
182  MOUT (output) INTEGER
183  If IJOB=1, the number of eigenvalues in the intervals.
184  If IJOB=2 or 3, the number of intervals output.
185  If IJOB=3, MOUT will equal MINP.
186 
187  NAB (input/output) INTEGER array, dimension (MMAX,2)
188  If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
189  If IJOB=2, then on input, NAB(i,j) should be set. It must
190  satisfy the condition:
191  N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
192  which means that in interval i only eigenvalues
193  NAB(i,1)+1,...,NAB(i,2) will be considered. Usually,
194  NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with
195  IJOB=1.
196  On output, NAB(i,j) will contain
197  max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
198  the input interval that the output interval
199  (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
200  the input values of NAB(k,1) and NAB(k,2).
201  If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
202  unless N(w) > NVAL(i) for all search points w , in which
203  case NAB(i,1) will not be modified, i.e., the output
204  value will be the same as the input value (modulo
205  reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
206  for all search points w , in which case NAB(i,2) will
207  not be modified. Normally, NAB should be set to some
208  distinctive value(s) before DLAEBZ is called.
209 
210  WORK (workspace) DOUBLE PRECISION array, dimension (MMAX)
211  Workspace.
212 
213  IWORK (workspace) INTEGER array, dimension (MMAX)
214  Workspace.
215 
216  INFO (output) INTEGER
217  = 0: All intervals converged.
218  = 1--MMAX: The last INFO intervals did not converge.
219  = MMAX+1: More than MMAX intervals were generated.
220 
221  Further Details
222  ===============
223 
224  This routine is intended to be called only by other LAPACK
225  routines, thus the interface is less user-friendly. It is intended
226  for two purposes:
227 
228  (a) finding eigenvalues. In this case, DLAEBZ should have one or
229  more initial intervals set up in AB, and DLAEBZ should be called
230  with IJOB=1. This sets up NAB, and also counts the eigenvalues.
231  Intervals with no eigenvalues would usually be thrown out at
232  this point. Also, if not all the eigenvalues in an interval i
233  are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
234  For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
235  eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX
236  no smaller than the value of MOUT returned by the call with
237  IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1
238  through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
239  tolerance specified by ABSTOL and RELTOL.
240 
241  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
242  In this case, start with a Gershgorin interval (a,b). Set up
243  AB to contain 2 search intervals, both initially (a,b). One
244  NVAL element should contain f-1 and the other should contain l
245  , while C should contain a and b, resp. NAB(i,1) should be -1
246  and NAB(i,2) should be N+1, to flag an error if the desired
247  interval does not lie in (a,b). DLAEBZ is then called with
248  IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals --
249  j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
250  if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
251  >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and
252  N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
253  w(l-r)=...=w(l+k) are handled similarly.
254 
255  =====================================================================
256 
257 
258  Check for Errors
259 
260  Parameter adjustments */
261  /* System generated locals */
262  integer nab_dim1, nab_offset, ab_dim1, ab_offset, i__1, i__2, i__3, i__4,
263  i__5, i__6;
264  Treal d__1, d__2, d__3, d__4;
265  /* Local variables */
266  integer itmp1, itmp2, j, kfnew, klnew, kf, ji, kl, jp, jit;
267  Treal tmp1, tmp2;
268 #define ab_ref(a_1,a_2) ab[(a_2)*ab_dim1 + a_1]
269 #define nab_ref(a_1,a_2) nab[(a_2)*nab_dim1 + a_1]
270 
271  nab_dim1 = *mmax;
272  nab_offset = 1 + nab_dim1 * 1;
273  nab -= nab_offset;
274  ab_dim1 = *mmax;
275  ab_offset = 1 + ab_dim1 * 1;
276  ab -= ab_offset;
277  --d__;
278  --e;
279  --e2;
280  --nval;
281  --c__;
282  --work;
283  --iwork;
284 
285  /* Function Body */
286  *info = 0;
287  if (*ijob < 1 || *ijob > 3) {
288  *info = -1;
289  return 0;
290  }
291 
292 /* Initialize NAB */
293 
294  if (*ijob == 1) {
295 
296 /* Compute the number of eigenvalues in the initial intervals. */
297 
298  *mout = 0;
299 /* DIR$ NOVECTOR */
300  i__1 = *minp;
301  for (ji = 1; ji <= i__1; ++ji) {
302  for (jp = 1; jp <= 2; ++jp) {
303  tmp1 = d__[1] - ab_ref(ji, jp);
304  if (absMACRO(tmp1) < *pivmin) {
305  tmp1 = -(*pivmin);
306  }
307  nab_ref(ji, jp) = 0;
308  if (tmp1 <= 0.) {
309  nab_ref(ji, jp) = 1;
310  }
311 
312  i__2 = *n;
313  for (j = 2; j <= i__2; ++j) {
314  tmp1 = d__[j] - e2[j - 1] / tmp1 - ab_ref(ji, jp);
315  if (absMACRO(tmp1) < *pivmin) {
316  tmp1 = -(*pivmin);
317  }
318  if (tmp1 <= 0.) {
319  nab_ref(ji, jp) = nab_ref(ji, jp) + 1;
320  }
321 /* L10: */
322  }
323 /* L20: */
324  }
325  *mout = *mout + nab_ref(ji, 2) - nab_ref(ji, 1);
326 /* L30: */
327  }
328  return 0;
329  }
330 
331 /* Initialize for loop
332 
333  KF and KL have the following meaning:
334  Intervals 1,...,KF-1 have converged.
335  Intervals KF,...,KL still need to be refined. */
336 
337  kf = 1;
338  kl = *minp;
339 
340 /* If IJOB=2, initialize C.
341  If IJOB=3, use the user-supplied starting point. */
342 
343  if (*ijob == 2) {
344  i__1 = *minp;
345  for (ji = 1; ji <= i__1; ++ji) {
346  c__[ji] = (ab_ref(ji, 1) + ab_ref(ji, 2)) * .5;
347 /* L40: */
348  }
349  }
350 
351 /* Iteration loop */
352 
353  i__1 = *nitmax;
354  for (jit = 1; jit <= i__1; ++jit) {
355 
356 /* Loop over intervals */
357 
358  if (kl - kf + 1 >= *nbmin && *nbmin > 0) {
359 
360 /* Begin of Parallel Version of the loop */
361 
362  i__2 = kl;
363  for (ji = kf; ji <= i__2; ++ji) {
364 
365 /* Compute N(c), the number of eigenvalues less than c */
366 
367  work[ji] = d__[1] - c__[ji];
368  iwork[ji] = 0;
369  if (work[ji] <= *pivmin) {
370  iwork[ji] = 1;
371 /* Computing MIN */
372  d__1 = work[ji], d__2 = -(*pivmin);
373  work[ji] = minMACRO(d__1,d__2);
374  }
375 
376  i__3 = *n;
377  for (j = 2; j <= i__3; ++j) {
378  work[ji] = d__[j] - e2[j - 1] / work[ji] - c__[ji];
379  if (work[ji] <= *pivmin) {
380  ++iwork[ji];
381 /* Computing MIN */
382  d__1 = work[ji], d__2 = -(*pivmin);
383  work[ji] = minMACRO(d__1,d__2);
384  }
385 /* L50: */
386  }
387 /* L60: */
388  }
389 
390  if (*ijob <= 2) {
391 
392 /* IJOB=2: Choose all intervals containing eigenvalues. */
393 
394  klnew = kl;
395  i__2 = kl;
396  for (ji = kf; ji <= i__2; ++ji) {
397 
398 /* Insure that N(w) is monotone
399 
400  Computing MIN
401  Computing MAX */
402  i__5 = nab_ref(ji, 1), i__6 = iwork[ji];
403  i__3 = nab_ref(ji, 2), i__4 = maxMACRO(i__5,i__6);
404  iwork[ji] = minMACRO(i__3,i__4);
405 
406 /* Update the Queue -- add intervals if both halves
407  contain eigenvalues. */
408 
409  if (iwork[ji] == nab_ref(ji, 2)) {
410 
411 /* No eigenvalue in the upper interval:
412  just use the lower interval. */
413 
414  ab_ref(ji, 2) = c__[ji];
415 
416  } else if (iwork[ji] == nab_ref(ji, 1)) {
417 
418 /* No eigenvalue in the lower interval:
419  just use the upper interval. */
420 
421  ab_ref(ji, 1) = c__[ji];
422  } else {
423  ++klnew;
424  if (klnew <= *mmax) {
425 
426 /* Eigenvalue in both intervals -- add upper to
427  queue. */
428 
429  ab_ref(klnew, 2) = ab_ref(ji, 2);
430  nab_ref(klnew, 2) = nab_ref(ji, 2);
431  ab_ref(klnew, 1) = c__[ji];
432  nab_ref(klnew, 1) = iwork[ji];
433  ab_ref(ji, 2) = c__[ji];
434  nab_ref(ji, 2) = iwork[ji];
435  } else {
436  *info = *mmax + 1;
437  }
438  }
439 /* L70: */
440  }
441  if (*info != 0) {
442  return 0;
443  }
444  kl = klnew;
445  } else {
446 
447 /* IJOB=3: Binary search. Keep only the interval containing
448  w s.t. N(w) = NVAL */
449 
450  i__2 = kl;
451  for (ji = kf; ji <= i__2; ++ji) {
452  if (iwork[ji] <= nval[ji]) {
453  ab_ref(ji, 1) = c__[ji];
454  nab_ref(ji, 1) = iwork[ji];
455  }
456  if (iwork[ji] >= nval[ji]) {
457  ab_ref(ji, 2) = c__[ji];
458  nab_ref(ji, 2) = iwork[ji];
459  }
460 /* L80: */
461  }
462  }
463 
464  } else {
465 
466 /* End of Parallel Version of the loop
467 
468  Begin of Serial Version of the loop */
469 
470  klnew = kl;
471  i__2 = kl;
472  for (ji = kf; ji <= i__2; ++ji) {
473 
474 /* Compute N(w), the number of eigenvalues less than w */
475 
476  tmp1 = c__[ji];
477  tmp2 = d__[1] - tmp1;
478  itmp1 = 0;
479  if (tmp2 <= *pivmin) {
480  itmp1 = 1;
481 /* Computing MIN */
482  d__1 = tmp2, d__2 = -(*pivmin);
483  tmp2 = minMACRO(d__1,d__2);
484  }
485 
486 /* A series of compiler directives to defeat vectorization
487  for the next loop
488 
489  $PL$ CMCHAR=' '
490  DIR$ NEXTSCALAR
491  $DIR SCALAR
492  DIR$ NEXT SCALAR
493  VD$L NOVECTOR
494  DEC$ NOVECTOR
495  VD$ NOVECTOR
496  VDIR NOVECTOR
497  VOCL LOOP,SCALAR
498  IBM PREFER SCALAR
499  $PL$ CMCHAR='*' */
500 
501  i__3 = *n;
502  for (j = 2; j <= i__3; ++j) {
503  tmp2 = d__[j] - e2[j - 1] / tmp2 - tmp1;
504  if (tmp2 <= *pivmin) {
505  ++itmp1;
506 /* Computing MIN */
507  d__1 = tmp2, d__2 = -(*pivmin);
508  tmp2 = minMACRO(d__1,d__2);
509  }
510 /* L90: */
511  }
512 
513  if (*ijob <= 2) {
514 
515 /* IJOB=2: Choose all intervals containing eigenvalues.
516 
517  Insure that N(w) is monotone
518 
519  Computing MIN
520  Computing MAX */
521  i__5 = nab_ref(ji, 1);
522  i__3 = nab_ref(ji, 2), i__4 = maxMACRO(i__5,itmp1);
523  itmp1 = minMACRO(i__3,i__4);
524 
525 /* Update the Queue -- add intervals if both halves
526  contain eigenvalues. */
527 
528  if (itmp1 == nab_ref(ji, 2)) {
529 
530 /* No eigenvalue in the upper interval:
531  just use the lower interval. */
532 
533  ab_ref(ji, 2) = tmp1;
534 
535  } else if (itmp1 == nab_ref(ji, 1)) {
536 
537 /* No eigenvalue in the lower interval:
538  just use the upper interval. */
539 
540  ab_ref(ji, 1) = tmp1;
541  } else if (klnew < *mmax) {
542 
543 /* Eigenvalue in both intervals -- add upper to queue. */
544 
545  ++klnew;
546  ab_ref(klnew, 2) = ab_ref(ji, 2);
547  nab_ref(klnew, 2) = nab_ref(ji, 2);
548  ab_ref(klnew, 1) = tmp1;
549  nab_ref(klnew, 1) = itmp1;
550  ab_ref(ji, 2) = tmp1;
551  nab_ref(ji, 2) = itmp1;
552  } else {
553  *info = *mmax + 1;
554  return 0;
555  }
556  } else {
557 
558 /* IJOB=3: Binary search. Keep only the interval
559  containing w s.t. N(w) = NVAL */
560 
561  if (itmp1 <= nval[ji]) {
562  ab_ref(ji, 1) = tmp1;
563  nab_ref(ji, 1) = itmp1;
564  }
565  if (itmp1 >= nval[ji]) {
566  ab_ref(ji, 2) = tmp1;
567  nab_ref(ji, 2) = itmp1;
568  }
569  }
570 /* L100: */
571  }
572  kl = klnew;
573 
574 /* End of Serial Version of the loop */
575 
576  }
577 
578 /* Check for convergence */
579 
580  kfnew = kf;
581  i__2 = kl;
582  for (ji = kf; ji <= i__2; ++ji) {
583  tmp1 = (d__1 = ab_ref(ji, 2) - ab_ref(ji, 1), absMACRO(d__1));
584 /* Computing MAX */
585  d__3 = (d__1 = ab_ref(ji, 2), absMACRO(d__1)), d__4 = (d__2 = ab_ref(
586  ji, 1), absMACRO(d__2));
587  tmp2 = maxMACRO(d__3,d__4);
588 /* Computing MAX */
589  d__1 = maxMACRO(*abstol,*pivmin), d__2 = *reltol * tmp2;
590  if (tmp1 < maxMACRO(d__1,d__2) || nab_ref(ji, 1) >= nab_ref(ji, 2)) {
591 
592 /* Converged -- Swap with position KFNEW,
593  then increment KFNEW */
594 
595  if (ji > kfnew) {
596  tmp1 = ab_ref(ji, 1);
597  tmp2 = ab_ref(ji, 2);
598  itmp1 = nab_ref(ji, 1);
599  itmp2 = nab_ref(ji, 2);
600  ab_ref(ji, 1) = ab_ref(kfnew, 1);
601  ab_ref(ji, 2) = ab_ref(kfnew, 2);
602  nab_ref(ji, 1) = nab_ref(kfnew, 1);
603  nab_ref(ji, 2) = nab_ref(kfnew, 2);
604  ab_ref(kfnew, 1) = tmp1;
605  ab_ref(kfnew, 2) = tmp2;
606  nab_ref(kfnew, 1) = itmp1;
607  nab_ref(kfnew, 2) = itmp2;
608  if (*ijob == 3) {
609  itmp1 = nval[ji];
610  nval[ji] = nval[kfnew];
611  nval[kfnew] = itmp1;
612  }
613  }
614  ++kfnew;
615  }
616 /* L110: */
617  }
618  kf = kfnew;
619 
620 /* Choose Midpoints */
621 
622  i__2 = kl;
623  for (ji = kf; ji <= i__2; ++ji) {
624  c__[ji] = (ab_ref(ji, 1) + ab_ref(ji, 2)) * .5;
625 /* L120: */
626  }
627 
628 /* If no more intervals to refine, quit. */
629 
630  if (kf > kl) {
631  goto L140;
632  }
633 /* L130: */
634  }
635 
636 /* Converged */
637 
638 L140:
639 /* Computing MAX */
640  i__1 = kl + 1 - kf;
641  *info = maxMACRO(i__1,0);
642  *mout = kl;
643 
644  return 0;
645 
646 /* End of DLAEBZ */
647 
648 } /* dlaebz_ */
649 
650 #undef nab_ref
651 #undef ab_ref
652 
653 
654 #endif
#define absMACRO(x)
Definition: template_blas_common.h:45
#define ab_ref(a_1, a_2)
int template_lapack_laebz(const integer *ijob, const integer *nitmax, const integer *n, const integer *mmax, const integer *minp, const integer *nbmin, const Treal *abstol, const Treal *reltol, const Treal *pivmin, const Treal *d__, const Treal *e, const Treal *e2, integer *nval, Treal *ab, Treal *c__, integer *mout, integer *nab, Treal *work, integer *iwork, integer *info)
Definition: template_lapack_laebz.h:40
int integer
Definition: template_blas_common.h:38
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
#define minMACRO(a, b)
Definition: template_blas_common.h:44
#define nab_ref(a_1, a_2)