454{
455 double* buff_2 = FLA_DOUBLE_PTR(
FLA_TWO );
456 double* buff_1 = FLA_DOUBLE_PTR(
FLA_ONE );
457 double* buff_0 = FLA_DOUBLE_PTR(
FLA_ZERO );
459
460 double first_elem;
461 double dot_product;
462 double beta, conj_beta;
463 double inv_tau11;
464 double minus_inv_tau11;
465 double minus_upsilon1, minus_conj_upsilon1;
466 double minus_psi1, minus_conj_psi1;
467 double minus_zeta1;
469
470
471 int b_alg = m_T;
472
473
474
475
476
477
478 double* buff_u = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
479 double* buff_y = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
480 double* buff_z = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
481 double* buff_v = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
482 double* buff_w = (
double* )
FLA_malloc( m_A *
sizeof( *buff_A ) );
483 int inc_u = 1;
484 int inc_y = 1;
485 int inc_z = 1;
486 int inc_v = 1;
487 int inc_w = 1;
488
489
490 first_elem = *buff_0;
491 minus_inv_tau11 = *buff_0;
492
493 for (
i = 0;
i < b_alg; ++
i )
494 {
495 double* A20 = buff_A + (0 )*cs_A + (
i+1)*rs_A;
496 double* alpha11 = buff_A + (
i )*cs_A + (
i )*rs_A;
497 double* a21 = buff_A + (
i )*cs_A + (
i+1)*rs_A;
498 double* A02 = buff_A + (
i+1)*cs_A + (0 )*rs_A;
499 double* a12t = buff_A + (
i+1)*cs_A + (
i )*rs_A;
500 double* A22 = buff_A + (
i+1)*cs_A + (
i+1)*rs_A;
501
502 double* t01 = buff_T + (
i )*cs_T + (0 )*rs_T;
503 double* tau11 = buff_T + (
i )*cs_T + (
i )*rs_T;
504
506 double* u2 = buff_u + (
i+1)*inc_u;
507
508 double* y0 = buff_y + (0 )*inc_y;
509 double*
psi1 = buff_y + (
i )*inc_y;
510 double* y2 = buff_y + (
i+1)*inc_y;
511
512 double*
zeta1 = buff_z + (
i )*inc_z;
513 double* z2 = buff_z + (
i+1)*inc_z;
514
515 double* v2 = buff_v + (
i+1)*inc_v;
516
517 double* w2 = buff_w + (
i+1)*inc_w;
518
519 double* a21_t = a21 + (0 )*cs_A + (0 )*rs_A;
520 double* a21_b = a21 + (0 )*cs_A + (1 )*rs_A;
521
522 int m_ahead = m_A -
i - 1;
523 int n_ahead = m_A -
i - 1;
526
527
528
529 if ( m_behind > 0 )
530 {
531
532
533
534 bl1_dmult3( buff_m1,
upsilon1, &minus_upsilon1 );
535 bl1_dcopyconj( &minus_upsilon1, &minus_conj_upsilon1 );
536
537
538
539
540 bl1_dmult3( buff_m1,
psi1, &minus_psi1 );
541 bl1_dcopyconj( &minus_psi1, &minus_conj_psi1 );
542
543
544
545 bl1_dmult3( buff_m1,
zeta1, &minus_zeta1 );
546
547
548
550 1,
551 &minus_upsilon1,
553 alpha11, 1 );
555 1,
556 &minus_zeta1,
558 alpha11, 1 );
559
560
561
563 m_ahead,
564 &minus_upsilon1,
565 y2, inc_y,
566 a12t, cs_A );
568 m_ahead,
569 &minus_zeta1,
570 u2, inc_u,
571 a12t, cs_A );
572
573
574
576 m_ahead,
577 &minus_conj_psi1,
578 u2, inc_u,
579 a21, rs_A );
581 m_ahead,
582 &minus_conj_upsilon1,
583 z2, inc_z,
584 a21, rs_A );
585 }
586
587 if ( m_ahead > 0 )
588 {
589
590
591
593 a21_t,
594 a21_b, rs_A,
595 tau11 );
596
597
598
599
600
601 bl1_ddiv3( buff_1, tau11, &inv_tau11 );
602 bl1_dneg2( &inv_tau11, &minus_inv_tau11 );
603
604
605
606 first_elem = *a21_t;
607 *a21_t = *buff_1;
608 }
609
610 if ( m_behind > 0 )
611 {
612
613
616 m_ahead,
617 n_ahead,
618 buff_m1,
619 u2, inc_u,
620 y2, inc_y,
621 A22, rs_A, cs_A );
624 m_ahead,
625 n_ahead,
626 buff_m1,
627 z2, inc_z,
628 u2, inc_u,
629 A22, rs_A, cs_A );
630 }
631
632 if ( m_ahead > 0 )
633 {
634
637 m_ahead,
638 n_ahead,
639 buff_1,
640 A22, rs_A, cs_A,
641 a21, rs_A,
642 buff_0,
643 v2, inc_v );
644
645
648 m_ahead,
649 n_ahead,
650 buff_1,
651 A22, rs_A, cs_A,
652 a21, rs_A,
653 buff_0,
654 w2, inc_w );
655
656
657
658
660 m_ahead,
661 a21, rs_A,
662 u2, inc_u );
664 m_ahead,
665 v2, inc_v,
666 y2, inc_y );
668 m_ahead,
669 w2, inc_w,
670 z2, inc_z );
671
672
673
674
676 m_ahead,
677 a21, rs_A,
678 z2, inc_z,
679 &beta );
680 bl1_dinvscals( buff_2, &beta );
681 bl1_dcopyconj( &beta, &conj_beta );
682
683
684
685
686 bl1_dscals( &minus_inv_tau11, &conj_beta );
688 m_ahead,
689 &conj_beta,
690 a21, rs_A,
691 y2, inc_y );
693 m_ahead,
694 &inv_tau11,
695 y2, inc_y );
696
697
698
699
700 bl1_dscals( &minus_inv_tau11, &beta );
702 m_ahead,
703 &beta,
704 a21, rs_A,
705 z2, inc_z );
707 m_ahead,
708 &inv_tau11,
709 z2, inc_z );
710
711
712
713
715 m_ahead,
716 a12t, cs_A,
717 a21, rs_A,
718 &dot_product );
719 bl1_dscals( &minus_inv_tau11, &dot_product );
721 m_ahead,
722 &dot_product,
723 a21, rs_A,
724 a12t, cs_A );
725
726
727
730 m_behind,
731 n_ahead,
732 buff_1,
733 A02, rs_A, cs_A,
734 a21, rs_A,
735 buff_0,
736 y0, inc_y );
739 m_behind,
740 n_ahead,
741 &minus_inv_tau11,
742 y0, inc_y,
743 a21, rs_A,
744 A02, rs_A, cs_A );
745
746
749 m_ahead,
750 n_behind,
751 buff_1,
752 A20, rs_A, cs_A,
753 a21, rs_A,
754 buff_0,
755 t01, rs_T );
756
757
758 *a21_t = first_elem;
759 }
760
761 if ( m_behind + 1 == b_alg && m_ahead > 0 )
762 {
763
764
767 m_ahead,
768 n_ahead,
769 buff_m1,
770 u2, inc_u,
771 y2, inc_y,
772 A22, rs_A, cs_A );
775 m_ahead,
776 n_ahead,
777 buff_m1,
778 z2, inc_z,
779 u2, inc_u,
780 A22, rs_A, cs_A );
781 }
782
783
784
785 }
786
787
788
789
790
791
797
798 return FLA_SUCCESS;
799}
FLA_Error FLA_Househ2_UT_l_opd(int m_x2, double *chi_1, double *x2, int inc_x2, double *tau)
Definition FLA_Househ2_UT.c:274
void bl1_daxpyv(conj1_t conj, int n, double *alpha, double *x, int incx, double *y, int incy)
Definition bl1_axpyv.c:21
void bl1_dcopyv(conj1_t conj, int m, double *x, int incx, double *y, int incy)
Definition bl1_copyv.c:42
void bl1_ddot(conj1_t conj, int n, double *x, int incx, double *y, int incy, double *rho)
Definition bl1_dot.c:26
void bl1_dgemv(trans1_t transa, conj1_t conjx, int m, int n, double *alpha, double *a, int a_rs, int a_cs, double *x, int incx, double *beta, double *y, int incy)
Definition bl1_gemv.c:69
void bl1_dger(conj1_t conjx, conj1_t conjy, int m, int n, double *alpha, double *x, int incx, double *y, int incy, double *a, int a_rs, int a_cs)
Definition bl1_ger.c:62
void bl1_dscalv(conj1_t conj, int n, double *alpha, double *x, int incx)
Definition bl1_scalv.c:24