<T>LAPACK 0.1.2
C++ Template Linear Algebra PACKage
Loading...
Searching...
No Matches
laed4.hpp
Go to the documentation of this file.
1
3//
4// Copyright (c) 2025 University of Colorado Denver. All rights reserved.
5//
6// This file is part of <T>LAPACK.
7// <T>LAPACK is free software: you can redistribute it and/or modify it under
8// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
9
10#ifndef TLAPACK_LAED4_HH
11#define TLAPACK_LAED4_HH
12
13//
14#include "tlapack/base/utils.hpp"
15
16//
17#include <tlapack/lapack/laed5.hpp>
18#include <tlapack/lapack/laed6.hpp>
19
20namespace tlapack {
21
93template <class d_t, class z_t, class delta_t, class real_t, class idx_t>
94int laed4(
95 idx_t n, idx_t i, d_t& d, z_t& z, delta_t& delta, real_t rho, real_t& dlam)
96{
97 int info = 0;
98 real_t psi, dpsi, phi, dphi, err, eta, a, b, c, w, del, tau, dltlb, dltub,
99 temp;
100
101 real_t maxIt = real_t(30);
102
103 if (n == 1) {
104 // Presumably, I = 1 upon entry
105 dlam = d[0] + rho * z[0] * z[0];
106 delta[0] = real_t(1.0);
107 return info;
108 }
109 else if (n == 2) {
110 laed5(i, d, z, delta, rho, dlam);
111 return info;
112 }
113
114 // Compute machine epsilon
115 real_t eps = ulp<real_t>();
116 // real_t eps = pow(2.0, -53);
117 real_t rhoinv = real_t(1.0) / rho;
118
119 // The Case if i = n
120 if (i == n - 1) {
121 // Initialize some basic variables
122 idx_t ii = n - 1;
123 idx_t niter = 1;
124
125 // Calculate Initial Guess
126 real_t midpt = rho / real_t(2.0);
127
128 // If ||Z||_2 is not one, then TEMP should be set to RHO * ||Z||_2^2 /
129 // TWO
130 for (idx_t j = 0; j < n; j++) {
131 delta[j] = (d[j] - d[i]) - midpt;
132 }
133
134 psi = real_t(0.0);
135 for (idx_t j = 0; j < n - 2; j++) {
136 psi += z[j] * z[j] / delta[j];
137 }
138
139 c = rhoinv + psi;
140 w = c + z[ii - 1] * z[ii - 1] / delta[ii - 1] +
141 z[n - 1] * z[n - 1] / delta[n - 1];
142
143 if (w <= 0) {
144 real_t temp = z[n - 2] * z[n - 2] / (d[n - 1] - d[n - 2] + rho) +
145 z[n - 1] * z[n - 1] / rho;
146
147 if (c <= temp) {
148 tau = rho;
149 }
150 else {
151 del = d[n - 1] - d[n - 2];
152 a = -c * del + z[n - 2] * z[n - 2] + z[n - 1] * z[n - 1];
153 b = z[n - 1] * z[n - 1] * del;
154
155 if (a < 0) {
156 tau = real_t(2.0) * b /
157 (sqrt(a * a + real_t(4.0) * b * c) - a);
158 }
159 else {
160 tau = (a + sqrt(a * a + real_t(4.0) * b * c)) /
161 (real_t(2.0) * c);
162 }
163 }
164
165 // It can be proved that D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <=
166 // D(N)+RHO
167
168 dltlb = midpt;
169 dltub = rho;
170 }
171 else {
172 del = d[n - 1] - d[n - 2];
173 a = -c * del + z[n - 2] * z[n - 2] + z[n - 1] * z[n - 1];
174 b = z[n - 1] * z[n - 1] * del;
175
176 if (a < 0)
177 tau = real_t(2.0) * b / (sqrt(a * a + real_t(4.0) * b * c) - a);
178 else
179 tau =
180 (a + sqrt(a * a + real_t(4.0) * b * c)) / (real_t(2.0) * c);
181
182 // It can be proved that* D(N) < D(N) + TAU < LAMBDA(N) < D(N) + RHO
183 // / 2 dltlb = 0.0;
184 dltlb = real_t(0.0);
185 dltub = midpt;
186 }
187
188 for (idx_t j = 0; j < n; j++) {
189 delta[j] = (d[j] - d[i]) - tau;
190 }
191
192 // Evaluate PSI and the derivative DPSI
193 dpsi = real_t(0.0);
194 phi = real_t(0.0);
195 dphi = real_t(0.0);
196 err = real_t(0.0);
197 psi = real_t(0.0);
198
199 for (idx_t j = 0; j < ii; j++) {
200 real_t temp = z[j] / delta[j];
201 psi += z[j] * temp;
202 dpsi += temp * temp;
203 err += psi;
204 }
205
206 err = abs(err);
207
208 // Evaluate PHI and the derivative DPHI
209 real_t temp = z[n - 1] / delta[n - 1];
210 phi = z[n - 1] * temp;
211 dphi = temp * temp;
212
213 err = real_t(8.0) * (-phi - psi) + err - phi + rhoinv +
214 abs(tau) * (dpsi + dphi);
215 w = rhoinv + phi + psi;
216
217 // Test for convergence
218 if (abs(w) <= eps * err) {
219 dlam = d[i] + tau;
220 return info;
221 }
222
223 if (w <= 0) {
224 dltlb = max(dltlb, tau);
225 }
226 else {
227 dltub = min(dltub, tau);
228 }
229
230 // Calculate the new step
231 niter++;
232 c = w - delta[n - 2] * dpsi - delta[n - 1] * dphi;
233 a = (delta[n - 2] + delta[n - 1]) * w -
234 delta[n - 2] * delta[n - 1] * (dpsi + dphi);
235 b = delta[n - 2] * delta[n - 1] * w;
236 if (c < 0) {
237 c = abs(c);
238 }
239 if (c == 0) {
240 // ETA = B/A
241 // ETA = RHO - TAU
242 // ETA = DLTUB - TAU
243 //
244 // Update proposed by Li, Ren-Cang:
245 eta = -w / (dpsi + dphi);
246 }
247 else if (a >= 0) {
248 eta = (a + sqrt(abs(a * a - real_t(4.0) * b * c))) /
249 (real_t(2.0) * c);
250 }
251 else {
252 eta =
253 real_t(2.0) * b / (a - sqrt(abs(a * a - real_t(4.0) * b * c)));
254 }
255
256 // Note, eta should be positive if w is negative, and
257 // eta should be negative otherwise. However,
258 // if for some reason caused by roundoff, eta*w > 0,
259 // we simply use one Newton step instead. This way
260 // will guarantee eta*w < 0.
261
262 if (w * eta > 0) {
263 eta = -w / (dpsi + dphi);
264 }
265
266 temp = tau + eta;
267 if (temp > dltub || temp < dltlb) {
268 if (w < 0) {
269 eta = (dltub - tau) / real_t(2.0);
270 }
271 else {
272 eta = (dltlb - tau) / real_t(2.0);
273 }
274 }
275
276 for (idx_t j = 0; j < n; j++) {
277 delta[j] -= eta;
278 }
279
280 tau += eta;
281
282 // Evaluate PSI and the derivative DPSI
283 dpsi = real_t(0.0);
284 psi = real_t(0.0);
285 err = real_t(0.0);
286 for (idx_t j = 0; j < ii; j++) {
287 temp = z[j] / delta[j];
288 psi += z[j] * temp;
289 dpsi += temp * temp;
290 err += psi;
291 }
292
293 err = abs(err);
294
295 // Evaluate PHI and the derivative DPHI
296 temp = z[n - 1] / delta[n - 1];
297 phi = z[n - 1] * temp;
298 dphi = temp * temp;
299 err = real_t(8.0) * (-phi - psi) + err - phi + rhoinv +
300 abs(tau) * (dpsi + dphi);
301
302 w = rhoinv + phi + psi;
303
304 // Main loop to update the values of the array DELTA
305
306 niter++;
307
308 while (niter < maxIt) {
309 // Test for convergence
310 if (abs(w) <= eps * err) {
311 dlam = d[i] + tau;
312 return info;
313 }
314
315 if (w <= 0) {
316 dltlb = max(dltlb, tau);
317 }
318 else {
319 dltub = min(dltub, tau);
320 }
321
322 // Calculate the new step
323 c = w - delta[n - 2] * dpsi - delta[n - 1] * dphi;
324 a = (delta[n - 2] + delta[n - 1]) * w -
325 delta[n - 2] * delta[n - 1] * (dpsi + dphi);
326 b = delta[n - 2] * delta[n - 1] * w;
327
328 if (a >= 0) {
329 eta = (a + sqrt(abs(a * a - real_t(4.0) * b * c))) /
330 (real_t(2.0) * c);
331 }
332 else {
333 eta = real_t(2.0) * b /
334 (a - sqrt(abs(a * a - real_t(4.0) * b * c)));
335 }
336
337 // Note, eta should be positive if w is negative, and eta should be
338 // negative otherwise. However, if for some reason caused by
339 // roundoff, eta*w > 0, we simply use one Newton step instead. This
340 // way will guarantee eta*w < 0.
341
342 if (w * eta > 0) {
343 eta = -w / (dpsi + dphi);
344 }
345 temp = tau + eta;
346 if (temp > dltub || temp < dltlb) {
347 if (w < 0) {
348 eta = (dltub - tau) / real_t(2.0);
349 }
350 else {
351 eta = (dltlb - tau) / real_t(2.0);
352 }
353 }
354 for (idx_t j = 0; j < n; j++) {
355 delta[j] -= eta;
356 }
357
358 tau = tau + eta;
359
360 // Evaluate PSI and the derivative DPSI
361 dpsi = real_t(0.0);
362 psi = real_t(0.0);
363 err = real_t(0.0);
364 for (idx_t j = 0; j < ii; j++) {
365 temp = z[j] / delta[j];
366 psi += z[j] * temp;
367 dpsi += temp * temp;
368 err = err + psi;
369 }
370 err = abs(err);
371
372 // Evaluate PHI and the derivative DPHI
373 temp = z[n - 1] / delta[n - 1];
374 phi = z[n - 1] * temp;
375 dphi = temp * temp;
376 err = real_t(8.0) * (-phi - psi) + err - phi + rhoinv +
377 abs(tau) * (dpsi + dphi);
378 w = rhoinv + phi + psi;
379
380 niter++;
381 }
382
383 // Return with INFO = 1, NITER = MAXIT and not converged
384 info = 1;
385 dlam = d[i] + tau;
386 return info;
387 }
388 else {
389 // The case for 0 ≤ i < n
390 idx_t niter = 0;
391 idx_t ip1 = i + 1;
392
393 // Calculate Inital Guess
394 del = d[ip1] - d[i];
395 real_t midpt = del / real_t(2.0);
396 for (idx_t j = 0; j < n; j++) {
397 delta[j] = (d[j] - d[i]) - midpt;
398 }
399
400 psi = real_t(0.0);
401 for (idx_t j = 0; j < i; j++) {
402 psi += z[j] * z[j] / delta[j];
403 }
404
405 phi = real_t(0.0);
406 for (idx_t j = n - 1; j >= i + 2; j--) {
407 phi += z[j] * z[j] / delta[j];
408 }
409
410 c = rhoinv + psi + phi;
411
412 w = c + z[i] * z[i] / delta[i] + z[ip1] * z[ip1] / delta[ip1];
413
414 bool orgati;
415
416 if (w > 0) {
417 // d(i)< the ith eigenvalue < (d(i)+d(i+1))/2
418 // We choose d(i) as origin.
419 orgati = true;
420 a = c * del + z[i] * z[i] + z[ip1] * z[ip1];
421
422 b = z[i] * z[i] * del;
423
424 if (a > 0) {
425 tau = real_t(2.0) * b /
426 (a + sqrt(abs(a * a - real_t(4.0) * b * c)));
427 }
428 else {
429 tau = (a - sqrt(abs(a * a - real_t(4.0) * b * c))) /
430 (real_t(2.0) * c);
431 }
432
433 dltlb = real_t(0.0);
434 dltub = midpt;
435 }
436 else {
437 // (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)
438 // We choose d(i+1) as origin.
439 orgati = false;
440 a = c * del - z[i] * z[i] - z[ip1] * z[ip1];
441 b = z[ip1] * z[ip1] * del;
442 if (a < 0) {
443 tau = real_t(2.0) * b /
444 (a - sqrt(abs(a * a + real_t(4.0) * b * c)));
445 }
446 else {
447 tau = -(a + sqrt(abs(a * a + real_t(4.0) * b * c))) /
448 (real_t(2.0) * c);
449 }
450
451 dltlb = -midpt;
452 dltub = real_t(0.0);
453 }
454
455 if (orgati) {
456 for (idx_t j = 0; j < n; j++) {
457 delta[j] = (d[j] - d[i]) - tau;
458 }
459 }
460 else {
461 for (idx_t j = 0; j < n; j++) {
462 delta[j] = (d[j] - d[ip1]) - tau;
463 }
464 }
465
466 idx_t ii;
467
468 if (orgati) {
469 ii = i;
470 }
471 else {
472 ii = i + 1;
473 }
474
475 idx_t iim1 = ii - 1;
476 idx_t iip1 = ii + 1;
477
478 // Evaluate PSI and the derivative DPSI
479 psi = real_t(0.0);
480 dpsi = real_t(0.0);
481 err = real_t(0.0);
482
483 for (idx_t j = 0; j + 1 <= ii; j++) {
484 temp = z[j] / delta[j];
485 psi = psi + z[j] * temp;
486 dpsi += temp * temp;
487 err += psi;
488 }
489 err = abs(err);
490
491 // Evaluate PHI and the derivative DPHI
492 phi = real_t(0.0);
493 dphi = real_t(0.0);
494
495 for (idx_t j = n - 1; j >= iip1; j--) {
496 temp = z[j] / delta[j];
497 phi += z[j] * temp;
498 dphi += temp * temp;
499 err += phi;
500 }
501
502 w = rhoinv + phi + psi;
503
504 // W is the value of the secular function with
505 // its ii-th element removed.
506
507 bool swtch3 = false;
508 if (orgati) {
509 if (w < 0) {
510 swtch3 = true;
511 }
512 }
513 else {
514 if (w > 0) {
515 swtch3 = true;
516 }
517 }
518
519 // if (ii == 1 || ii == n) {
520 if (ii == 0 || ii == n - 1) {
521 swtch3 = false;
522 }
523
524 temp = z[ii] / delta[ii];
525 real_t dw = dpsi + dphi + temp * temp;
526 temp = z[ii] * temp;
527 w += temp;
528 err = real_t(8.0) * (phi - psi) + err + real_t(2.0) * rhoinv +
529 real_t(3.0) * abs(temp) + abs(tau) * dw;
530
531 // Test for Convergence
532 if (abs(w) <= eps * err) {
533 if (orgati) {
534 dlam = d[i] + tau;
535 }
536 else {
537 dlam = d[ip1] + tau;
538 }
539 return info;
540 }
541
542 if (w <= 0) {
543 dltlb = max(dltlb, tau);
544 }
545 else {
546 dltub = min(dltub, tau);
547 }
548
549 // Calculate the new step
550 niter++;
551
552 real_t zz[3];
553
554 if (!swtch3) {
555 if (orgati) {
556 c = w - delta[ip1] * dw -
557 (d[i] - d[ip1]) * ((z[i] / delta[i]) * (z[i] / delta[i]));
558 }
559 else {
560 c = w - delta[i] * dw -
561 (d[ip1] - d[i]) * ((z[ip1]) / delta[ip1]) *
562 ((z[ip1]) / delta[ip1]);
563 }
564 a = (delta[i] + delta[ip1]) * w - delta[i] * delta[ip1] * dw;
565 b = delta[i] * delta[ip1] * w;
566 if (c == 0) {
567 if (a == 0) {
568 if (orgati) {
569 a = z[i] * z[i] +
570 delta[ip1] * delta[ip1] * (dpsi + dphi);
571 }
572 else {
573 a = z[ip1] * z[ip1] +
574 delta[i] * delta[i] * (dpsi + dphi);
575 }
576 }
577 eta = b / a;
578 }
579 else if (a <= 0) {
580 eta = (a - sqrt(abs(a * a - real_t(4.0) * b * c))) /
581 (real_t(2.0) * c);
582 }
583 else {
584 eta = real_t(2.0) * b /
585 (a + sqrt(abs(a * a - real_t(4.0) * b * c)));
586 }
587 }
588 else {
589 // Interpolation using THREE most relevant poles
590 temp = rhoinv + psi + phi;
591 if (orgati) {
592 real_t temp1 = z[iim1] / delta[iim1];
593 temp1 = temp1 * temp1;
594 c = temp - delta[iip1] * (dpsi + dphi) -
595 (d[iim1] - d[iip1]) * temp1;
596 zz[0] = z[iim1] * z[iim1];
597 zz[2] = delta[iip1] * delta[iip1] * ((dpsi - temp1) + dphi);
598 }
599 else {
600 real_t temp1 = z[iip1] / delta[iip1];
601 temp1 = temp1 * temp1;
602 c = temp - delta[iim1] * (dpsi + dphi) -
603 (d[iip1] - d[iim1]) * temp1;
604 zz[0] = delta[iim1] * delta[iim1] * (dpsi + (dphi - temp1));
605 zz[2] = z[iip1] * z[iip1];
606 }
607 zz[1] = z[ii] * z[ii];
608
609 std::vector<real_t> sub(delta.begin() + iim1, delta.end());
610 info = laed6(niter, orgati, c, sub, zz, w, eta);
611
612 if (info == 0) {
613 return info;
614 }
615 }
616
617 // Note, eta should be positive if w is negative, and
618 // eta should be negative otherwise. However,
619 // if for some reason caused by roundoff, eta*w > 0,
620 // we simply use one Newton step instead. This way
621 // will guarantee eta*w < 0.
622
623 if (w * eta >= 0) {
624 eta = -w / dw;
625 }
626 temp = tau + eta;
627 if (temp > dltub || temp < dltlb) {
628 if (w < 0) {
629 eta = (dltub - tau) / real_t(2.0);
630 }
631 else {
632 eta = (dltlb - tau) / real_t(2.0);
633 }
634 }
635
636 real_t prew = w;
637
638 for (idx_t j = 0; j < n; j++) {
639 delta[j] -= eta;
640 }
641
642 // Evaluate PSI and the derivative DPSI
643 psi = real_t(0.0);
644 dpsi = real_t(0.0);
645 err = real_t(0.0);
646
647 for (idx_t j = 0; j + 1 <= ii; j++) {
648 temp = z[j] / delta[j];
649 psi += z[j] * temp;
650 dpsi += temp * temp;
651 err += psi;
652 }
653 err = abs(err);
654
655 // Evaluate PHI and the derivative DPHI
656 phi = real_t(0.0);
657 dphi = real_t(0.0);
658 for (idx_t j = n - 1; j >= iip1; j--) {
659 temp = z[j] / delta[j];
660 phi += z[j] * temp;
661 dphi += temp * temp;
662 err += phi;
663 }
664
665 temp = z[ii] / delta[ii];
666 dw = dpsi + dphi + temp * temp;
667 temp = z[ii] * temp;
668 w = rhoinv + phi + psi + temp;
669 err = real_t(8.0) * (phi - psi) + err + real_t(2.0) * rhoinv +
670 real_t(3.0) * abs(temp) + abs(tau + eta) * dw;
671
672 bool swtch = false;
673 if (orgati) {
674 if (-w > abs(prew) / real_t(10.0)) {
675 swtch = true;
676 }
677 }
678 else {
679 if (w > abs(prew) / real_t(10.0)) {
680 swtch = true;
681 }
682 }
683
684 tau = tau + eta;
685
686 // Main loop to update the values of the array DELTA
687
688 idx_t iter = niter + 1;
689
690 while (iter < maxIt) {
691 // Test for convergence
692 if (abs(w) <= eps * err) {
693 if (orgati) {
694 dlam = d[i] + tau;
695 }
696 else {
697 dlam = d[ip1] + tau;
698 }
699 return info;
700 }
701
702 if (w <= 0) {
703 dltlb = max(dltlb, tau);
704 }
705 else {
706 dltub = min(dltub, tau);
707 }
708
709 // Calculate the new step
710 if (!swtch3) {
711 if (!swtch) {
712 if (orgati) {
713 c = w - delta[ip1] * dw -
714 (d[i] - d[ip1]) * (z[i] / delta[i]) *
715 (z[i] / delta[i]);
716 }
717 else {
718 c = w - delta[i] * dw -
719 (d[ip1] - d[i]) * (z[ip1] / delta[ip1]) *
720 (z[ip1] / delta[ip1]);
721 }
722 }
723 else {
724 temp = z[ii] / delta[ii];
725 if (orgati) {
726 dpsi += temp * temp;
727 }
728 else {
729 dphi += temp * temp;
730 }
731 c = w - delta[i] * dpsi - delta[ip1] * dphi;
732 }
733
734 a = (delta[i] + delta[ip1]) * w - delta[i] * delta[ip1] * dw;
735 b = delta[i] * delta[ip1] * w;
736
737 if (c == 0) {
738 if (a == 0) {
739 if (!swtch) {
740 if (orgati) {
741 a = z[i] * z[i] +
742 delta[ip1] * delta[ip1] * (dpsi + dphi);
743 }
744 else {
745 a = z[ip1] * z[ip1] +
746 delta[i] * delta[i] * (dpsi + dphi);
747 }
748 }
749 else {
750 a = delta[i] * delta[i] * dpsi +
751 delta[ip1] * delta[ip1] * dphi;
752 }
753 }
754
755 eta = b / a;
756 }
757 else if (a < 0) {
758 eta = (a - sqrt(abs(a * a - real_t(4.0) * b * c))) /
759 (real_t(2.0) * c);
760 }
761 else {
762 eta = real_t(2.0) * b /
763 (a + sqrt(abs(a * a - real_t(4.0) * b * c)));
764 }
765 }
766 else {
767 // Interpolation using THREE most relevant poles
768 temp = rhoinv + psi + phi;
769 if (swtch) {
770 c = temp - delta[iim1] * dpsi - delta[iip1] * dphi;
771 zz[0] = delta[iim1] * delta[iim1] * dpsi;
772 zz[2] = delta[iip1] * delta[iip1] * dphi;
773 }
774 else {
775 if (orgati) {
776 real_t temp1 = z[iim1] / delta[iim1];
777 temp1 = temp1 * temp1;
778 c = temp - delta[iip1] * (dpsi + dphi) -
779 (d[iim1] - d[iip1]) * temp1;
780 zz[0] = z[iim1] * z[iim1];
781 zz[2] =
782 delta[iip1] * delta[iip1] * ((dpsi - temp1) + dphi);
783 }
784 else {
785 real_t temp1 = z[iip1] / delta[iip1];
786 temp1 = temp1 * temp1;
787 c = temp - delta[iim1] * (dpsi + dphi) -
788 (d[iip1] - d[iim1]) * temp1;
789 zz[0] =
790 delta[iim1] * delta[iim1] * (dpsi + (dphi - temp1));
791 zz[2] = z[iip1] * z[iip1];
792 }
793 }
794
795 std::vector<real_t> sub(delta.begin() + iim1, delta.end());
796 info = laed6(niter, orgati, c, sub, zz, w, eta);
797
798 if (info == 0) {
799 return info;
800 }
801 }
802
803 // Note, eta should be positive if w is negative, and
804 // eta should be negative otherwise. However,
805 // if for some reason caused by roundoff, eta*w > 0,
806 // we simply use one Newton step instead. This way
807 // will guarantee eta*w < 0.
808
809 if (w * eta >= 0) {
810 eta = -w / dw;
811 }
812 temp = tau + eta;
813 if (temp > dltub || temp < dltlb) {
814 if (w < 0) {
815 eta = (dltub - tau) / real_t(2.0);
816 }
817 else {
818 eta = (dltlb - tau) / real_t(2.0);
819 }
820 }
821
822 for (idx_t j = 0; j < n; j++) {
823 delta[j] = delta[j] - eta;
824 }
825
826 tau += eta;
827 prew = w;
828
829 // Evaluate PSI and the derivative DPSI
830 psi = real_t(0.0);
831 dpsi = real_t(0.0);
832 err = real_t(0.0);
833 for (idx_t j = 0; j + 1 <= ii; j++) {
834 temp = z[j] / delta[j];
835 psi += z[j] * temp;
836 dpsi += temp * temp;
837 err += psi;
838 }
839
840 err = abs(err);
841
842 // Evaluate PHI and the derivative DPHI
843 phi = real_t(0.0);
844 dphi = real_t(0.0);
845 for (idx_t j = n - 1; j >= iip1; j--) {
846 temp = z[j] / delta[j];
847 phi += z[j] * temp;
848 dphi += temp * temp;
849 err += phi;
850 }
851
852 temp = z[ii] / delta[ii];
853 dw = dpsi + dphi + temp * temp;
854 temp = z[ii] * temp;
855 w = rhoinv + phi + psi + temp;
856 err = real_t(8.0) * (phi - psi) + err + real_t(2.0) * rhoinv +
857 real_t(3.0) * abs(temp) + abs(tau) * dw;
858
859 if (w * prew > 0 && abs(w) > abs(prew) / real_t(10.0)) {
860 swtch = !swtch;
861 }
862
863 iter++;
864 }
865
866 info = 1;
867 if (orgati) {
868 dlam = d[i] + tau;
869 }
870 else {
871 dlam = d[ip1] + tau;
872 }
873 }
874
875 return info;
876}
877} // namespace tlapack
878
879#endif // TLAPACK_LAED4_HH
tlapack
Sort the numbers in D in increasing order (if ID = 'I') or in decreasing order (if ID = 'D' ).
Definition arrayTraits.hpp:15
tlapack::real_type
typename traits::real_type_traits< Types..., int >::type real_type
The common real type of the list of types.
Definition scalar_type_traits.hpp:113
tlapack::laed6
int laed6(idx_t kniter, bool &orgati, real_t rho, d_t &d, z_t &z, real_t &finit, real_t &tau)
LAED6 used by STEDC.
Definition laed6.hpp:78
tlapack::laed5
void laed5(idx_t i, d_t &d, z_t &z, delta_t &delta, real_t rho, real_t &dlam)
LAED5 used by STEDC.
Definition laed5.hpp:61
tlapack::laed4
int laed4(idx_t n, idx_t i, d_t &d, z_t &z, delta_t &delta, real_t rho, real_t &dlam)
LAED4 used by STEDC.
Definition laed4.hpp:94