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