14 #include "../PBpblas.h"
15 #include "../PBtools.h"
16 #include "../PBblacs.h"
17 #include "../PBblas.h"
21 char * DIAG,
int N,
char * A,
int IA,
int JA,
22 int * DESCA,
char * XC,
int INCXC,
char * XR,
25 void PB_Cptrsv(
TYPE, FBCAST, UPLO, TRANS, DIAG, N, A, IA, JA, DESCA,
26 XC, INCXC, XR, INCXR )
30 char * DIAG, * TRANS, * UPLO;
31 int FBCAST, IA, INCXC, INCXR, JA, N;
229 char btop, * negone, * one, * zero;
230 int Acol, Aii, Aimb1, Ainb1, Ais1Col, Ais1Row, AisColRep,
231 AisRowRep, Ajj, Alcol, Ald, Alrow, Amb, Anpprev, Anb, Anp,
232 Anq, Arow, Asrc, ChangeRoc=0, bcst, ctxt, ione=1, k=0, kb,
233 kbprev=0, kbsize, mb1, mycol, myrow, n1, n1last, n1p,
234 n1pprev=0, nb1, nlast, npcol, nprow, rocprev, size, tmp1,
248 char * Aprev = NULL, * Xd = NULL, * Xdprev = NULL,
249 * Xprev = NULL, * work = NULL;
262 PB_Cinfog2l( IA, JA, DESCA, nprow, npcol, myrow, mycol, &Aii, &Ajj, &Arow,
268 Amb = DESCA[
MB_]; Anb = DESCA[
NB_]; Ald = DESCA[
LLD_ ];
270 Anp =
PB_Cnumroc( N, 0, Aimb1, Amb, myrow, Arow, nprow );
271 Ais1Row = !(
PB_Cspan( N, 0, Aimb1, Amb, Arow, nprow ) );
273 Anq =
PB_Cnumroc( N, 0, Ainb1, Anb, mycol, Acol, npcol );
274 Ais1Col = !(
PB_Cspan( N, 0, Ainb1, Anb, Acol, npcol ) );
278 if( Ais1Row && Ais1Col )
288 TYPE->size ), &Ald, XC, &INCXC );
289 TYPE->Fcopy( &Anp, XC, &INCXC, XR, &INCXR );
291 if( ( Arow >= 0 ) && FBCAST )
295 TYPE->Cgebs2d( ctxt,
COLUMN, &btop, 1, Anq, XR, INCXR );
297 TYPE->Cgebr2d( ctxt,
COLUMN, &btop, 1, Anq, XR, INCXR, Arow,
310 TYPE->size ), &Ald, XR, &INCXR );
311 TYPE->Fcopy( &Anq, XR, &INCXR, XC, &INCXC );
313 if( Acol >= 0 && FBCAST )
317 TYPE->Cgebs2d( ctxt,
ROW, &btop, Anp, 1, XC, Anp );
319 TYPE->Cgebr2d( ctxt,
ROW, &btop, Anp, 1, XC, Anp, myrow,
330 negone =
TYPE->negone; one =
TYPE->one; zero =
TYPE->zero;
331 axpy =
TYPE->Faxpy; copy =
TYPE->Fcopy; set =
TYPE->Fset;
332 gemv =
TYPE->Fgemv; trsv =
TYPE->Ftrsv;
333 send =
TYPE->Cgesd2d; recv =
TYPE->Cgerv2d;
334 bsend =
TYPE->Cgebs2d; brecv =
TYPE->Cgebr2d;
336 if( ( Anp > 0 ) && ( Anq > 0 ) ) A =
Mptr( A, Aii, Ajj, Ald, size );
340 if( ( Anq <= 0 ) || ( Ais1Row && ( ( Arow >= 0 ) && !( FBCAST ) &&
341 ( myrow != Arow ) ) ) )
return;
343 bcst = ( ( !Ais1Row ) || ( Ais1Row && ( Arow >= 0 ) && FBCAST ) );
344 AisRowRep = ( ( Arow < 0 ) || ( nprow == 1 ) );
351 nlast = ( npcol - 1 ) * Anb;
352 n1 =
MAX( nlast, Anb );
354 n1last = n1 - Anb +
MAX( Ainb1, Anb );
357 Alrow =
PB_Cindxg2p( tmp1, Aimb1, Amb, Arow, Arow, nprow );
358 Alcol =
PB_Cindxg2p( tmp1, Ainb1, Anb, Acol, Acol, npcol );
359 rocprev = Alcol; Anpprev = Anp; Xprev = XC; Xdprev = XR;
360 Aprev = A =
Mptr( A, 0, Anq, Ald, size );
363 tmp1 = N - ( kb =
MIN( mb1, nb1 ) );
364 n1 = ( ( Ais1Col || ( N - nb1 < nlast ) ) ? n1last : n1 );
365 tmp2 = n1 + nb1 - kb; tmp1 -= ( tmp2 =
MIN( tmp1, tmp2 ) );
367 n1p =
PB_Cnumroc( tmp2,
MAX( 0, tmp1 ), Aimb1, Amb, myrow, Asrc,
373 if( Ais1Col || ( mycol == Alcol ) )
377 Xd =
Mptr( XR, 0, Anq, INCXR, size );
379 if( ( Arow < 0 ) || ( myrow == Alrow ) ) { Anp -= kb; }
385 if( ( Ais1Col || ( mycol == rocprev ) ) && ( kbprev > 0 ) )
387 tmp1 = ( Anpprev - n1pprev ) * size;
388 gemv(
C2F_CHAR( TRANS ), &n1pprev, &kbprev, negone,
389 Aprev+tmp1, &Ald, Xdprev, &INCXR, one, Xprev+tmp1,
395 if( !( Ais1Col ) && ChangeRoc )
397 if( mycol == rocprev )
399 send( ctxt, n1pprev, 1, Xprev+(Anpprev-n1pprev)*size,
400 n1pprev, myrow, Alcol );
402 else if( mycol == Alcol )
404 recv( ctxt, n1pprev, 1, work, n1pprev, myrow, rocprev );
405 axpy( &n1pprev, one, work, &ione,
Mptr( Xprev,
406 Anpprev-n1pprev, 0, INCXC, size ), &INCXC );
413 if( Ais1Col || ( mycol == Alcol ) )
415 if( AisRowRep || ( myrow == Alrow ) )
418 &kb,
Mptr( A, Anp, 0, Ald, size ), &Ald,
Mptr( XC, Anp,
419 0, INCXC, size ), &INCXC );
420 copy( &kb,
Mptr( XC, Anp, 0, INCXC, size ), &INCXC,
Mptr( XR,
421 0, Anq, INCXR, size ), &INCXR );
426 bsend( ctxt,
COLUMN, &btop, 1, kb,
Mptr( XR, 0, Anq, INCXR,
429 brecv( ctxt,
COLUMN, &btop, 1, kb,
Mptr( XR, 0, Anq, INCXR,
430 size ), INCXR, Alrow, mycol );
435 if( !( Ais1Col ) && ( AisRowRep || ( myrow == Alrow ) ) )
436 set( &kb, zero,
Mptr( XC, Anp, 0, INCXC, size ), &ione );
441 if( ( Ais1Col || ( mycol == rocprev ) ) && ( kbprev > 0 ) &&
442 ( ( tmp1 = Anpprev - n1pprev ) > 0 ) )
443 gemv(
C2F_CHAR( TRANS ), &tmp1, &kbprev, negone, Aprev, &Ald,
444 Xdprev, &INCXR, one, Xprev, &INCXC );
448 if( Ais1Col || ( mycol == Alcol ) ) { Xdprev = Xd; Aprev = A; }
449 if( AisRowRep || ( myrow == Alrow ) ) { Anpprev -= kb; }
460 if( !( Ais1Row ) && ( Alrow >= 0 ) )
462 mb1 = ( N > Aimb1 ? Amb : Aimb1 );
466 ChangeRoc = ( nb1 == 0 );
470 if( !( Ais1Col ) && ( Alcol >= 0 ) )
472 nb1 = ( N > Ainb1 ? Anb : Ainb1 );
474 tmp1 = N - ( kb =
MIN( mb1, nb1 ) );
475 n1 = ( ( Ais1Col || ( N - nb1 < nlast ) ) ? n1last : n1 );
476 tmp2 = n1 + nb1 - kb; tmp1 -= ( tmp2 =
MIN( tmp1, tmp2 ) );
477 n1p =
PB_Cnumroc( tmp2,
MAX( 0, tmp1 ), Aimb1, Amb, myrow, Asrc,
486 n1 = (
MAX( npcol, 2 ) - 1 ) * Anb;
488 Aprev = A; Xprev = XC; Xdprev = XR; Anpprev = Anp;
489 mb1 = Aimb1; nb1 = Ainb1; rocprev = Acol;
490 tmp1 = N - ( kb =
MIN( mb1, nb1 ) ); tmp2 = n1 + nb1 - kb;
492 n1p =
PB_Cnumroc(
MIN( tmp1, tmp2 ), kb, Aimb1, Amb, myrow, Asrc,
502 if( ( Ais1Col || ( mycol == rocprev ) ) && ( kbprev > 0 ) )
503 gemv(
C2F_CHAR( TRANS ), &n1pprev, &kbprev, negone, Aprev,
504 &Ald, Xdprev, &INCXR, one, Xprev, &INCXC );
508 if( !( Ais1Col ) && ChangeRoc )
510 if( mycol == rocprev )
512 send( ctxt, n1pprev, 1, Xprev, n1pprev, myrow, Acol );
514 else if( mycol == Acol )
516 recv( ctxt, n1pprev, 1, work, n1pprev, myrow, rocprev );
517 axpy( &n1pprev, one, work, &ione, Xprev, &INCXC );
524 if( Ais1Col || ( mycol == Acol ) )
526 if( AisRowRep || ( myrow == Arow ) )
529 &kb, A, &Ald, XC, &INCXC );
530 copy( &kb, XC, &INCXC, XR, &INCXR );
535 bsend( ctxt,
COLUMN, &btop, 1, kb, XR, INCXR );
537 brecv( ctxt,
COLUMN, &btop, 1, kb, XR, INCXR, Arow,
543 if( !( Ais1Col ) && ( AisRowRep || ( myrow == Arow ) ) )
544 set( &kb, zero, XC, &INCXC );
549 if( ( Ais1Col || ( mycol == rocprev ) ) && ( kbprev > 0 ) )
551 if( ( tmp1 = Anpprev - n1pprev ) > 0 )
553 tmp2 = n1pprev * size;
554 gemv(
C2F_CHAR( TRANS ), &tmp1, &kbprev, negone, Aprev+tmp2,
555 &Ald, Xdprev, &INCXR, one, Xprev+tmp2, &INCXC );
557 Aprev += Ald * kbprev * size;
562 if( Ais1Col || ( mycol == Acol ) )
563 { A += Ald*kbsize; Xdprev = Xd = XR; XR += INCXR*kbsize; }
564 if( AisRowRep || ( myrow == Arow ) )
566 Xprev = ( XC += kbsize );
569 Anpprev = ( Anp -= kb );
580 if( !( Ais1Row ) && ( Arow >= 0 ) )
586 ChangeRoc = ( nb1 == 0 );
590 if( !( Ais1Col ) && ( Acol >= 0 ) )
594 tmp1 = N - ( kb =
MIN( mb1, nb1 ) ); tmp2 = n1 + nb1 - kb;
595 n1p =
PB_Cnumroc(
MIN( tmp2, tmp1 ), k+kb, Aimb1, Amb, myrow, Asrc,
602 if( ( Anp <= 0 ) || ( Ais1Col && ( ( Acol >= 0 ) && !( FBCAST ) &&
603 ( mycol != Acol ) ) ) )
return;
605 bcst = ( ( !Ais1Col ) || ( Ais1Col && ( Acol >= 0 ) && FBCAST ) );
606 AisColRep = ( ( Acol < 0 ) || ( npcol == 1 ) );
613 n1 = (
MAX( nprow, 2 ) - 1 ) * Amb;
615 Aprev = A; Xprev = XR; Xdprev = XC; Anpprev = Anq;
616 mb1 = Aimb1; nb1 = Ainb1; rocprev = Arow;
617 tmp1 = N - ( kb =
MIN( mb1, nb1 ) ); tmp2 = n1 + mb1 - kb;
619 n1p =
PB_Cnumroc(
MIN( tmp1, tmp2 ), kb, Ainb1, Anb, mycol, Asrc,
629 if( ( Ais1Row || ( myrow == rocprev ) ) && ( kbprev > 0 ) )
630 gemv(
C2F_CHAR( TRANS ), &kbprev, &n1pprev, negone, Aprev,
631 &Ald, Xdprev, &INCXC, one, Xprev, &INCXR );
635 if( !( Ais1Row ) && ChangeRoc )
637 if( myrow == rocprev )
639 send( ctxt, 1, n1pprev, Xprev, INCXR, Arow, mycol );
641 else if( myrow == Arow )
643 recv( ctxt, 1, n1pprev, work, 1, rocprev, mycol );
644 axpy( &n1pprev, one, work, &ione, Xprev, &INCXR );
651 if( Ais1Row || ( myrow == Arow ) )
653 if( AisColRep || ( mycol == Acol ) )
656 &kb, A, &Ald, XR, &INCXR );
657 copy( &kb, XR, &INCXR, XC, &INCXC );
662 bsend( ctxt,
ROW, &btop, kb, 1, XC, kb );
664 brecv( ctxt,
ROW, &btop, kb, 1, XC, kb, myrow, Acol );
669 if( !( Ais1Row ) && ( AisColRep || ( mycol == Acol ) ) )
670 set( &kb, zero, XR, &INCXR );
675 if( ( Ais1Row || ( myrow == rocprev ) ) && ( kbprev > 0 ) )
677 if( ( tmp1 = Anpprev - n1pprev ) > 0 )
679 tmp2 = n1pprev * size;
680 gemv(
C2F_CHAR( TRANS ), &kbprev, &tmp1, negone,
681 Aprev+Ald*tmp2, &Ald, Xdprev, &INCXC, one,
682 Xprev+INCXR*tmp2, &INCXR );
684 Aprev += kbprev * size;
689 if( Ais1Row || ( myrow == Arow ) )
690 { A += kbsize; Xdprev = Xd = XC; XC += kbsize; }
691 if( AisColRep || ( mycol == Acol ) )
693 Xprev = ( XR += INCXR * kbsize );
695 Anpprev = ( Anq -= kb );
696 Aprev += Ald * kbsize;
707 if( !( Ais1Col ) && ( Acol >= 0 ) )
713 ChangeRoc = ( mb1 == 0 );
717 if( !( Ais1Row ) && ( Arow >= 0 ) )
721 tmp1 = N - ( kb =
MIN( mb1, nb1 ) ); tmp2 = n1 + mb1 - kb;
722 n1p =
PB_Cnumroc(
MIN( tmp2, tmp1 ), k+kb, Ainb1, Anb, mycol, Asrc,
731 nlast = ( nprow - 1 ) * Amb;
732 n1 =
MAX( nlast, Amb );
734 n1last = n1 - Amb +
MAX( Aimb1, Amb );
737 Alrow =
PB_Cindxg2p( tmp1, Aimb1, Amb, Arow, Arow, nprow );
738 Alcol =
PB_Cindxg2p( tmp1, Ainb1, Anb, Acol, Acol, npcol );
739 rocprev = Alrow; Anpprev = Anq; Xprev = XR; Xdprev = XC;
740 Aprev = A =
Mptr( A, Anp, 0, Ald, size );
743 tmp1 = N - ( kb =
MIN( mb1, nb1 ) );
744 n1 = ( ( Ais1Row || ( N - mb1 < nlast ) ) ? n1last : n1 );
745 tmp2 = n1 + mb1 - kb; tmp1 -= ( tmp2 =
MIN( tmp1, tmp2 ) );
747 n1p =
PB_Cnumroc( tmp2,
MAX( 0, tmp1 ), Ainb1, Anb, mycol, Asrc,
753 if( Ais1Row || ( myrow == Alrow ) )
757 Xd =
Mptr( XC, Anp, 0, INCXC, size );
759 if( ( Acol < 0 ) || ( mycol == Alcol ) ) { Anq -= kb; }
765 if( ( Ais1Row || ( myrow == rocprev ) ) && ( kbprev > 0 ) )
767 tmp1 = ( Anpprev - n1pprev ) * size;
768 gemv(
C2F_CHAR( TRANS ), &kbprev, &n1pprev, negone,
769 Aprev+Ald*tmp1, &Ald, Xdprev, &INCXC, one,
770 Xprev+INCXR*tmp1, &INCXR );
775 if( !( Ais1Row ) && ChangeRoc )
777 if( myrow == rocprev )
779 send( ctxt, 1, n1pprev,
Mptr( Xprev, 0, Anpprev-n1pprev,
780 INCXR, size ), INCXR, Alrow, mycol );
782 else if( myrow == Alrow )
784 recv( ctxt, 1, n1pprev, work, 1, rocprev, mycol );
785 axpy( &n1pprev, one, work, &ione,
Mptr( Xprev, 0,
786 Anpprev-n1pprev, INCXR, size ), &INCXR );
793 if( Ais1Row || ( myrow == Alrow ) )
795 if( AisColRep || ( mycol == Alcol ) )
798 &kb,
Mptr( A, 0, Anq, Ald, size ), &Ald,
Mptr( XR, 0,
799 Anq, INCXR, size ), &INCXR );
800 copy( &kb,
Mptr( XR, 0, Anq, INCXR, size ), &INCXR,
Mptr( XC,
801 0, Anp, INCXC, size ), &INCXC );
806 bsend( ctxt,
ROW, &btop, kb, 1,
Mptr( XC, 0, Anp, INCXC,
809 brecv( ctxt,
ROW, &btop, kb, 1,
Mptr( XC, 0, Anp, INCXC,
810 size ), kb, myrow, Alcol );
815 if( !( Ais1Row ) && ( AisColRep || ( mycol == Alcol ) ) )
816 set( &kb, zero,
Mptr( XR, 0, Anq, INCXR, size ), &INCXR );
821 if( ( Ais1Row || ( myrow == rocprev ) ) && ( kbprev > 0 ) &&
822 ( ( tmp1 = Anpprev - n1pprev ) > 0 ) )
823 gemv(
C2F_CHAR( TRANS ), &kbprev, &tmp1, negone, Aprev, &Ald,
824 Xdprev, &INCXC, one, Xprev, &INCXR );
828 if( Ais1Row || ( myrow == Alrow ) ) { Xdprev = Xd; Aprev = A; }
829 if( AisColRep || ( mycol == Alcol ) ) { Anpprev -= kb; }
840 if( !( Ais1Col ) && ( Alcol >= 0 ) )
842 nb1 = ( N > Ainb1 ? Anb : Ainb1 );
846 ChangeRoc = ( mb1 == 0 );
850 if( !( Ais1Row ) && ( Alrow >= 0 ) )
852 mb1 = ( N > Aimb1 ? Amb : Aimb1 );
854 tmp1 = N - ( kb =
MIN( mb1, nb1 ) );
855 n1 = ( ( Ais1Row || ( N - mb1 < nlast ) ) ? n1last : n1 );
856 tmp2 = n1 + mb1 - kb; tmp1 -= ( tmp2 =
MIN( tmp1, tmp2 ) );
857 n1p =
PB_Cnumroc( tmp2,
MAX( 0, tmp1 ), Ainb1, Anb, mycol, Asrc,
862 if( work ) free( work );