From 4adaeaf5d547d7c092102aa9f84278b47fa1b4c6 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 16 Mar 2006 17:36:34 +0000 Subject: [PATCH] Fixed bug in matrix regeneration stuff. --- src/serial/dp/dcocr.f | 359 +++++++++++++++++++------------------- src/serial/psb_cest.f90 | 8 +- src/serial/psb_dcoins.f90 | 10 +- src/serial/psb_dcsdp.f90 | 91 ++++++++-- src/tools/psb_dspasb.f90 | 49 +++--- 5 files changed, 292 insertions(+), 225 deletions(-) diff --git a/src/serial/dp/dcocr.f b/src/serial/dp/dcocr.f index 941967a8..7c484cd2 100644 --- a/src/serial/dp/dcocr.f +++ b/src/serial/dp/dcocr.f @@ -41,7 +41,7 @@ C C C .. Scalar Arguments .. INTEGER LARN, LAUX, LAUX2, LIAN1, LIAN2, M, - + N, IUPDUP, IERROR + + N, IUPDUP, IERROR CHARACTER TRANS,UNITD C .. Array Arguments .. DOUBLE PRECISION AR(*), ARN(*), D(*) @@ -74,7 +74,7 @@ C CALL FCPSB_ERRACTIONSAVE(ERR_ACT) CHECK_FLAG=IBITS(INFO(PSB_UPD_),1,2) -c$$$ write(0,*) 'DCOCR FLAG ',info(upd_),check_flag +c$$$ write(0,*) 'DCOCR FLAG ',info(psb_upd_),check_flag IF (TRANS.EQ.'N') THEN SCALE = (UNITD.EQ.'L') ! meaningless @@ -112,8 +112,8 @@ C C Error handling C IF(IERROR.NE.0) THEN - CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) - GOTO 9999 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 END IF IF (DESCRA(1:1).EQ.'G') THEN @@ -126,12 +126,15 @@ c$$$ do k=1,nnz c$$$ write(*,*) k,ia(k),ja(k),ar(k) c$$$ enddo c$$$ endif +c$$$ write(0,*) 'DCOCR Sizes ',lian2,((m+1)+nnz+psb_ireg_flgs_+1), +c$$$ + m+1,nnz,psb_ireg_flgs_, +c$$$ + laux,2*(2+nnz) if ((lian2.ge.((m+1)+nnz+psb_ireg_flgs_+1)) + .and.(laux.ge.2*(2+nnz))) then C C Prepare for smart regeneration c - + ipx = nnz+3 do i=1, nnz aux(ipx+i-1) = i @@ -145,7 +148,10 @@ c ian2(ip1+psb_nnz_) = 0 ian2(ip1+psb_ichk_) = nnz+check_flag -c$$$ write(0,*)'DCOCR m,ip1,ip2,nnz',m, +c$$$ write(0,*)'DCOCR Check: ',ip2,ian2(ip1+psb_iflag_), +c$$$ + ian2(ip1+psb_nnzt_), ian2(ip1+psb_nnz_), +c$$$ + ian2(ip1+psb_ichk_) + c$$$ + ip1,ip2,nnz,ian2(ip1+nnzt_) if (debug) write(0,*) 'Build check :',ian2(ip1+psb_nnzt_) @@ -153,12 +159,12 @@ C .... Order with key IA ... CALL MRGSRT(NNZ,IA,AUX,IRET) IF (IRET.EQ.0) CALL REORDVN3(NNZ,AR,IA,JA,AUX(IPX),AUX) if (debug) then - do i=1, nnz-1 - if (ia(i).gt.ia(i+1)) then - write(0,*) 'Sorting error:',i,ia(i),ia(i+1) - endif - enddo - write(0,*) 'nnz :',m,nnz,ia(nnz),ja(nnz) + do i=1, nnz-1 + if (ia(i).gt.ia(i+1)) then + write(0,*) 'Sorting error:',i,ia(i),ia(i+1) + endif + enddo + write(0,*) 'nnz :',m,nnz,ia(nnz),ja(nnz) endif C .... Order with key IA2N ... @@ -217,9 +223,9 @@ C ... Insert other element of row ... ELSE IF (CHECK_FLAG.EQ.1) THEN C ... Error, there are duplicated elements ... - IERROR = 130 - CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) - GOTO 9999 + IERROR = 130 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 ELSE IF (CHECK_FLAG.EQ.2) THEN C ... Insert only the last duplicated element ... ARN(ELEM_CSR-1) = AR(ELEM) @@ -227,7 +233,7 @@ C ... Insert only the last duplicated element ... ELSE IF (CHECK_FLAG.EQ.3) THEN C ... Sum the duplicated element ... ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM) - ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1 + ian2(ip2+aux(ipx+elem-1)-1) = elem_csr-1 END IF ENDIF ELEM = ELEM + 1 @@ -235,6 +241,7 @@ C ... Sum the duplicated element ... IAN2(ROW+1) = ELEM_CSR ENDDO ELSE +c$$$ write(0,*) 'Going for ELSE !!!?' C .... Order with key IA ... CALL MRGSRT(NNZ,IA,AUX,IRET) IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) @@ -291,19 +298,19 @@ C ... Insert other element of row ... ELSE IF (CHECK_FLAG.EQ.1) THEN C ... Error, there are duplicated elements ... - IERROR = 130 - CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) - GOTO 9999 + IERROR = 130 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 ELSE IF (CHECK_FLAG.EQ.2) THEN C ... Insert only the last duplicated element ... ARN(ELEM_CSR-1) = AR(ELEM) if (debug) write(0,*) 'Duplicated overwrite', - + elem_csr-1,elem + + elem_csr-1,elem ELSE IF (CHECK_FLAG.EQ.3) THEN C ... Sum the duplicated element ... ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM) if (debug) write(0,*) 'Duplicated add', - + elem_csr-1,elem + + elem_csr-1,elem END IF ENDIF ELEM = ELEM + 1 @@ -313,11 +320,11 @@ C ... Sum the duplicated element ... ENDIF if (debug) write(0,*)'Done Rebuild CSR', - + ian2(m+1),ia(elem) + + ian2(m+1),ia(elem) if (debug) then - do i=ian2(m+1), nnz - write(0,*) 'Overflow check :',ia(i),ja(i),ar(i) - enddo + do i=ian2(m+1), nnz + write(0,*) 'Overflow check :',ia(i),ja(i),ar(i) + enddo endif ELSE IF (DESCRA(1:1).EQ.'S' .AND. DESCRA(2:2).EQ.'U') THEN @@ -329,197 +336,197 @@ C ... Sum the duplicated element ... ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'U') THEN C .... Order with key IA ... - CALL MRGSRT(NNZ,IA,AUX,IRET) - IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) + CALL MRGSRT(NNZ,IA,AUX,IRET) + IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) C .... Order with key IA2N ... - I = 1 - J = I + I = 1 + J = I c$$$ DO WHILE (I.LE.NNZ) c$$$ DO WHILE ((IA(J).EQ.IA(I)).AND. c$$$ + (J.LE.NNZ)) - DO - if (I>NNZ) exit - DO - if (j>nnz) exit - if (ia(j) /= ia(i)) exit - J = J+1 - ENDDO - NZL = J - I - CALL MRGSRT(NZL,JA(I),AUX,IRET) - IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX) - I = J + DO + if (I>NNZ) exit + DO + if (j>nnz) exit + if (ia(j) /= ia(i)) exit + J = J+1 ENDDO + NZL = J - I + CALL MRGSRT(NZL,JA(I),AUX,IRET) + IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX) + I = J + ENDDO C ... Construct CSR Representation... - ELEM = 1 - ELEM_CSR = 1 + ELEM = 1 + ELEM_CSR = 1 C ... Insert first element ... - DO ROW = 1, IA(1) - IAN2(ROW) = 1 - ENDDO - if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr - IF(JA(ELEM).GT.IA(ELEM)) THEN - IAN1(ELEM_CSR) = JA(ELEM) - ARN(ELEM_CSR) = AR(ELEM) - ELEM_CSR = ELEM_CSR+1 - ENDIF + DO ROW = 1, IA(1) + IAN2(ROW) = 1 + ENDDO + if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr + IF(JA(ELEM).GT.IA(ELEM)) THEN + IAN1(ELEM_CSR) = JA(ELEM) + ARN(ELEM_CSR) = AR(ELEM) + ELEM_CSR = ELEM_CSR+1 + ENDIF - ELEM = ELEM+1 + ELEM = ELEM+1 C ... Insert remaining element ... - DO ROW = IA(1), M + DO ROW = IA(1), M c$$$ if (debug) write(*,*)'CSR Loop:',row,m,elem_csr c$$$ DO WHILE ((IA(ELEM).EQ.ROW).AND.(ELEM.LE.NNZ)) - DO - if (elem > nnz) exit - if (ia(elem) /= row) exit - IF (IA(ELEM).NE.IA(ELEM-1)) THEN + DO + if (elem > nnz) exit + if (ia(elem) /= row) exit + IF (IA(ELEM).NE.IA(ELEM-1)) THEN C ... Insert first element of a row ... - IF(JA(ELEM).GT.IA(ELEM)) THEN - IAN1(ELEM_CSR) = JA(ELEM) - ARN(ELEM_CSR) = AR(ELEM) - ELEM_CSR = ELEM_CSR+1 - ENDIF - ELSE IF (JA(ELEM).NE.JA(ELEM-1)) THEN + IF(JA(ELEM).GT.IA(ELEM)) THEN + IAN1(ELEM_CSR) = JA(ELEM) + ARN(ELEM_CSR) = AR(ELEM) + ELEM_CSR = ELEM_CSR+1 + ENDIF + ELSE IF (JA(ELEM).NE.JA(ELEM-1)) THEN C ... Insert other element of row ... - IF(JA(ELEM).GT.IA(ELEM)) THEN - IAN1(ELEM_CSR) = JA(ELEM) - ARN(ELEM_CSR) = AR(ELEM) - ELEM_CSR = ELEM_CSR+1 - ENDIF - ELSE - IF (CHECK_FLAG.EQ.1) THEN + IF(JA(ELEM).GT.IA(ELEM)) THEN + IAN1(ELEM_CSR) = JA(ELEM) + ARN(ELEM_CSR) = AR(ELEM) + ELEM_CSR = ELEM_CSR+1 + ENDIF + ELSE + IF (CHECK_FLAG.EQ.1) THEN C ... Error, there are duplicated elements ... - IERROR = 130 - CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) - GOTO 9999 - ELSE IF (CHECK_FLAG.EQ.2) THEN + IERROR = 130 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 + ELSE IF (CHECK_FLAG.EQ.2) THEN C ... Insert only the last duplicated element ... - IF(JA(ELEM).GT.IA(ELEM)) THEN - ARN(ELEM_CSR-1) = AR(ELEM) - ENDIF - if (debug) write(0,*) 'Duplicated overwrite', - + elem_csr-1,elem - ELSE IF (CHECK_FLAG.EQ.3) THEN + IF(JA(ELEM).GT.IA(ELEM)) THEN + ARN(ELEM_CSR-1) = AR(ELEM) + ENDIF + if (debug) write(0,*) 'Duplicated overwrite', + + elem_csr-1,elem + ELSE IF (CHECK_FLAG.EQ.3) THEN C ... Sum the duplicated element ... - IF(JA(ELEM).GT.IA(ELEM)) THEN - ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM) - ENDIF - if (debug) write(0,*) 'Duplicated add', - + elem_csr-1,elem - END IF + IF(JA(ELEM).GT.IA(ELEM)) THEN + ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM) ENDIF - ELEM = ELEM + 1 - ENDDO - IAN2(ROW+1) = ELEM_CSR + if (debug) write(0,*) 'Duplicated add', + + elem_csr-1,elem + END IF + ENDIF + ELEM = ELEM + 1 ENDDO + IAN2(ROW+1) = ELEM_CSR + ENDDO - - if (debug) write(0,*)'Done Rebuild CSR', - + ian2(m+1),ia(elem) - if (debug) then - do i=ian2(m+1), nnz - write(0,*) 'Overflow check :',ia(i),ja(i),ar(i) - enddo - endif - - - - ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') THEN + + if (debug) write(0,*)'Done Rebuild CSR', + + ian2(m+1),ia(elem) + if (debug) then + do i=ian2(m+1), nnz + write(0,*) 'Overflow check :',ia(i),ja(i),ar(i) + enddo + endif + + + + ELSE IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') THEN C .... Order with key IA ... - CALL MRGSRT(NNZ,IA,AUX,IRET) - IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) + CALL MRGSRT(NNZ,IA,AUX,IRET) + IF (IRET.EQ.0) CALL REORDVN(NNZ,AR,IA,JA,AUX) C .... Order with key IA2N ... - I = 1 - J = I + I = 1 + J = I c$$$ DO WHILE (I.LE.NNZ) c$$$ DO WHILE ((IA(J).EQ.IA(I)).AND. c$$$ + (J.LE.NNZ)) - DO - if (I>NNZ) exit - DO - if (j>nnz) exit - if (ia(j) /= ia(i)) exit - J = J+1 - ENDDO - NZL = J - I - CALL MRGSRT(NZL,JA(I),AUX,IRET) - IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX) - I = J + DO + if (I>NNZ) exit + DO + if (j>nnz) exit + if (ia(j) /= ia(i)) exit + J = J+1 ENDDO + NZL = J - I + CALL MRGSRT(NZL,JA(I),AUX,IRET) + IF (IRET.EQ.0) CALL REORDVN(NZL,AR(I),IA(I),JA(I),AUX) + I = J + ENDDO C ... Construct CSR Representation... - ELEM = 1 - ELEM_CSR = 1 + ELEM = 1 + ELEM_CSR = 1 C ... Insert first element ... - DO ROW = 1, IA(1) - IAN2(ROW) = 1 - ENDDO - if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr - IF(JA(ELEM).LT.IA(ELEM)) THEN - IAN1(ELEM_CSR) = JA(ELEM) - ARN(ELEM_CSR) = AR(ELEM) - ELEM_CSR = ELEM_CSR+1 - ENDIF - ELEM = ELEM+1 - + DO ROW = 1, IA(1) + IAN2(ROW) = 1 + ENDDO + if (debug) write(0,*)'Rebuild CSR',ia(1),elem_csr + IF(JA(ELEM).LT.IA(ELEM)) THEN + IAN1(ELEM_CSR) = JA(ELEM) + ARN(ELEM_CSR) = AR(ELEM) + ELEM_CSR = ELEM_CSR+1 + ENDIF + ELEM = ELEM+1 + C ... Insert remaining element ... - DO ROW = IA(1), M + DO ROW = IA(1), M c$$$ if (debug) write(*,*)'CSR Loop:',row,m,elem_csr c$$$ DO WHILE ((IA(ELEM).EQ.ROW).AND.(ELEM.LE.NNZ)) - DO - if (elem > nnz) exit - if (ia(elem) /= row) exit - IF (IA(ELEM).NE.IA(ELEM-1)) THEN + DO + if (elem > nnz) exit + if (ia(elem) /= row) exit + IF (IA(ELEM).NE.IA(ELEM-1)) THEN C ... Insert first element of a row ... - IF(JA(ELEM).LT.IA(ELEM)) THEN - IAN1(ELEM_CSR) = JA(ELEM) - ARN(ELEM_CSR) = AR(ELEM) - ELEM_CSR = ELEM_CSR+1 - ENDIF - ELSE IF (JA(ELEM).NE.JA(ELEM-1)) THEN + IF(JA(ELEM).LT.IA(ELEM)) THEN + IAN1(ELEM_CSR) = JA(ELEM) + ARN(ELEM_CSR) = AR(ELEM) + ELEM_CSR = ELEM_CSR+1 + ENDIF + ELSE IF (JA(ELEM).NE.JA(ELEM-1)) THEN C ... Insert other element of row ... - IF(JA(ELEM).LT.IA(ELEM)) THEN - IAN1(ELEM_CSR) = JA(ELEM) - ARN(ELEM_CSR) = AR(ELEM) - ELEM_CSR = ELEM_CSR+1 - ENDIF - ELSE - IF (CHECK_FLAG.EQ.1) THEN + IF(JA(ELEM).LT.IA(ELEM)) THEN + IAN1(ELEM_CSR) = JA(ELEM) + ARN(ELEM_CSR) = AR(ELEM) + ELEM_CSR = ELEM_CSR+1 + ENDIF + ELSE + IF (CHECK_FLAG.EQ.1) THEN C ... Error, there are duplicated elements ... - IERROR = 130 - CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) - GOTO 9999 - ELSE IF (CHECK_FLAG.EQ.2) THEN + IERROR = 130 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 + ELSE IF (CHECK_FLAG.EQ.2) THEN C ... Insert only the last duplicated element ... - IF(JA(ELEM).LT.IA(ELEM)) THEN - ARN(ELEM_CSR-1) = AR(ELEM) - ENDIF - if (debug) write(0,*) 'Duplicated overwrite', - + elem_csr-1,elem - ELSE IF (CHECK_FLAG.EQ.3) THEN + IF(JA(ELEM).LT.IA(ELEM)) THEN + ARN(ELEM_CSR-1) = AR(ELEM) + ENDIF + if (debug) write(0,*) 'Duplicated overwrite', + + elem_csr-1,elem + ELSE IF (CHECK_FLAG.EQ.3) THEN C ... Sum the duplicated element ... - IF(JA(ELEM).LT.IA(ELEM)) THEN - ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM) - ENDIF - if (debug) write(0,*) 'Duplicated add', - + elem_csr-1,elem - END IF - ENDIF - ELEM = ELEM + 1 - ENDDO - IAN2(ROW+1) = ELEM_CSR + IF(JA(ELEM).LT.IA(ELEM)) THEN + ARN(ELEM_CSR-1) = ARN(ELEM_CSR-1) + AR(ELEM) + ENDIF + if (debug) write(0,*) 'Duplicated add', + + elem_csr-1,elem + END IF + ENDIF + ELEM = ELEM + 1 ENDDO + IAN2(ROW+1) = ELEM_CSR + ENDDO if (debug) write(0,*)'Done Rebuild CSR', - + ian2(m+1),ia(elem) + + ian2(m+1),ia(elem) if (debug) then - do i=ian2(m+1), nnz - write(0,*) 'Overflow check :',ia(i),ja(i),ar(i) - enddo + do i=ian2(m+1), nnz + write(0,*) 'Overflow check :',ia(i),ja(i),ar(i) + enddo endif @@ -529,9 +536,9 @@ C C C TO DO C - IERROR = 3021 - CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) - GOTO 9999 + IERROR = 3021 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 END IF INFON(1)=ELEM_CSR-1 @@ -543,8 +550,8 @@ C CALL FCPSB_ERRACTIONRESTORE(ERR_ACT) IF ( ERR_ACT .NE. 0 ) THEN - CALL FCPSB_SERROR() - RETURN + CALL FCPSB_SERROR() + RETURN ENDIF RETURN diff --git a/src/serial/psb_cest.f90 b/src/serial/psb_cest.f90 index 9977520a..3a1f3175 100644 --- a/src/serial/psb_cest.f90 +++ b/src/serial/psb_cest.f90 @@ -28,14 +28,14 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_cest(afmt, nnz, lia1, lia2, lar, up, info) +subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info) use psb_error_mod use psb_const_mod implicit none ! .. scalar arguments .. - integer :: nnz, lia1, lia2, lar, info + integer :: m,n,nnz, lia1, lia2, lar, info character :: up ! .. array arguments.. character(len=5) :: afmt @@ -52,7 +52,7 @@ subroutine psb_cest(afmt, nnz, lia1, lia2, lar, up, info) if ((up.eq.'y').or.(up.eq.'Y')) then if (afmt.eq.'JAD') then lia1 = 2*(nnz + nnz/5) +1000 - lia2 = 2*(nnz + nnz/5) +1000 + lia2 = 2*(nnz + nnz/5) +1000 +m lar = nnz + nnz/5 else if (afmt.eq.'COO') then lia1 = nnz @@ -60,7 +60,7 @@ subroutine psb_cest(afmt, nnz, lia1, lia2, lar, up, info) lar = nnz else if(afmt.eq.'CSR') then lia1 = nnz - lia2 = 2*nnz + 1000 + lia2 = 2*nnz + 1000 + m + 1 lar = nnz else info = 136 diff --git a/src/serial/psb_dcoins.f90 b/src/serial/psb_dcoins.f90 index 5a32d0cf..d4d2f687 100644 --- a/src/serial/psb_dcoins.f90 +++ b/src/serial/psb_dcoins.f90 @@ -149,11 +149,11 @@ subroutine psb_dcoins(nz,ia,ja,val,a,gtl,imin,imax,jmin,jmax,info) call psb_errpush(info,name,a_err=ch_err) goto 9999 endif -!!$ if (debug) then -!!$ if ((nza - a%ia2(ip1+nnz_)) /= nz) then -!!$ write(0,*) 'PSB_COINS: update discarded items ' -!!$ end if -!!$ end if + if (debug) then + if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then + write(0,*) 'PSB_COINS: update discarded items ' + end if + end if a%ia2(ip1+psb_nnz_) = nza else diff --git a/src/serial/psb_dcsdp.f90 b/src/serial/psb_dcsdp.f90 index ca4c64df..e4bbb522 100644 --- a/src/serial/psb_dcsdp.f90 +++ b/src/serial/psb_dcsdp.f90 @@ -61,13 +61,31 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd) real(kind(1.d0)), allocatable :: work(:) type(psb_dspmat_type) :: temp_a Integer :: nzr, ntry, ifc_,ierror, ia1_size,& - & ia2_size, aspk_size + & ia2_size, aspk_size,size_req,n_row,n_col,iup integer :: ip1, ip2, nnz, iflag, ichk, nnzt,& & ipc, i, count, err_act, ierrv(5) - character :: check_,trans_,unitd_ + character :: check_,trans_,unitd_, up Integer, Parameter :: maxtry=8 logical, parameter :: debug=.false. character(len=20) :: name, ch_err + interface psb_cest + subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info) + integer, intent(in) :: m,n,nnz + integer, intent(out) :: lia1, lia2, lar, info + character, intent(inout) :: afmt*5 + character, intent(in) :: up + end subroutine psb_cest + end interface + + interface psb_spinfo + subroutine psb_dspinfo(ireq,a,ires,info,iaux) + use psb_spmat_type + type(psb_dspmat_type), intent(in) :: a + integer, intent(in) :: ireq + integer, intent(out) :: ires, info + integer, intent(in), optional :: iaux + end subroutine psb_dspinfo + end interface name='psb_dcsdp' info = 0 @@ -133,6 +151,45 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd) ! ...matrix conversion... b%m=a%m b%k=a%k + call psb_spinfo(psb_nztotreq_,a,size_req,info) + if (debug) write(0,*) 'DCSDP : size_req 1:',size_req + !! PULL IUP FROM INFOA FIELD + iup = iand(b%infoa(psb_upd_),4) + if (iup > 0) then + up = 'Y' + else + up = 'N' + endif + n_row=b%m + n_col=b%k + call psb_cest(b%fida, n_row,n_col,size_req,& + & ia1_size, ia2_size, aspk_size, up,info) + +!!$ write(0,*) 'ESTIMATE : ',ia1_size,ia2_size,aspk_Size,iup + if (info /= no_err) then + info=4010 + ch_err='psb_cest' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + if ((size(b%aspk) < aspk_size) .or.& + &(size(b%ia1) < ia1_size) .or.& + &(size(b%ia2) < ia2_size) .or.& + &(size(b%pl) < b%m) .or.& + &(size(b%pr) < b%k )) then + call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) + endif + if (info /= no_err) then + info=4010 + ch_err='psb_sp_reall' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + endif + + b%pl(:) = 0 + b%pr(:) = 0 + + select case (a%fida(1:3)) case ('CSR') @@ -141,11 +198,11 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd) case ('CSR') - - ia1_size=a%infoa(psb_nnz_) - ia2_size=a%m+1 - aspk_size=a%infoa(psb_nnz_) - call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) +!!$ +!!$ ia1_size=a%infoa(psb_nnz_) +!!$ ia2_size=a%m+1 +!!$ aspk_size=a%infoa(psb_nnz_) +!!$ call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) call dcrcr(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& & a%ia1, a%ia2, a%infoa, b%pl, b%descra, b%aspk, b%ia1,& @@ -164,10 +221,10 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd) !...converting to JAD !...output matrix may not be big enough - ia1_size=a%infoa(psb_nnz_) - ia2_size=a%m+1 - aspk_size=a%infoa(psb_nnz_) - call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) +!!$ ia1_size=a%infoa(psb_nnz_) +!!$ ia2_size=a%m+1 +!!$ aspk_size=a%infoa(psb_nnz_) +!!$ call psb_sp_reall(b,ia1_size,ia2_size,aspk_size,info) do call dcrjd(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& @@ -207,8 +264,8 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd) case ('COO') - aspk_size=max(size(a%aspk),a%ia2(a%m+1)) - call psb_sp_reall(b,aspk_size,info) +!!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1)) +!!$ call psb_sp_reall(b,aspk_size,info) !!$ write(0,*) 'From DCSDP90:',b%fida,size(b%aspk),info call dcrco(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& & a%ia1, a%ia2, a%infoa, b%pl, b%descra, b%aspk, b%ia1,& @@ -228,8 +285,8 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd) case ('CSR') - aspk_size=max(size(a%aspk),a%ia2(a%m+1)) - call psb_sp_reall(b,aspk_size,info) +!!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1)) +!!$ call psb_sp_reall(b,aspk_size,info) call dcocr(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& & a%ia2, a%ia1, a%infoa, b%pl, b%descra, b%aspk, b%ia1,& & b%ia2, b%infoa, b%pr, size(b%aspk), size(b%ia1),& @@ -301,8 +358,8 @@ subroutine psb_dcsdp(a, b,info,ifc,check,trans,unitd) case ('COO') - aspk_size=max(size(a%aspk),a%ia2(a%m+1)) - call psb_sp_reall(b,aspk_size,info) +!!$ aspk_size=max(size(a%aspk),a%ia2(a%m+1)) +!!$ call psb_sp_reall(b,aspk_size,info) call dcoco(trans_, a%m, a%k, unitd_, d, a%descra, a%aspk,& & a%ia1, a%ia2, a%infoa, b%pl, b%descra, b%aspk, b%ia1,& & b%ia2, b%infoa, b%pr, size(b%aspk), size(b%ia1),& diff --git a/src/tools/psb_dspasb.f90 b/src/tools/psb_dspasb.f90 index a83c35f5..71fb0c88 100644 --- a/src/tools/psb_dspasb.f90 +++ b/src/tools/psb_dspasb.f90 @@ -54,10 +54,11 @@ subroutine psb_dspasb(a,desc_a, info, afmt, up, dup) implicit none interface psb_cest - subroutine psb_cest(afmt, nnz, lia1, lia2, lar, up, info) - integer, intent(in) :: nnz + subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, up, info) + integer, intent(in) :: m,n,nnz integer, intent(out) :: lia1, lia2, lar, info - character, intent(in) :: afmt*5, up + character, intent(inout) :: afmt*5 + character, intent(in) :: up end subroutine psb_cest end interface @@ -179,26 +180,28 @@ subroutine psb_dspasb(a,desc_a, info, afmt, up, dup) ! work area requested must be fixed to ! No of Grid'd processes and NNZ+2 ! - size_req = max(a%infoa(psb_nnz_),1)+3 - if (debug) write(0,*) 'DCSDP : size_req 1:',size_req - call psb_cest(a%fida, size_req, ia1_size, ia2_size, aspk_size, iup,info) - if (info /= no_err) then - info=4010 - ch_err='psb_cest' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - call psb_sp_reall(a,ia1_size,ia2_size,aspk_size,info) - if (info /= no_err) then - info=4010 - ch_err='psb_sp_reall' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - endif - - a%pl(:) = 0 - a%pr(:) = 0 +!!$ size_req = max(a%infoa(psb_nnz_),1)+3 +!!$ if (debug) write(0,*) 'DCSDP : size_req 1:',size_req +!!$ call psb_cest(a%fida, n_row,n_col,size_req,& +!!$ & ia1_size, ia2_size, aspk_size, iup,info) +!!$ write(0,*) 'ESTIMATE : ',ia1_size,ia2_size,aspk_Size,iup +!!$ if (info /= no_err) then +!!$ info=4010 +!!$ ch_err='psb_cest' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ endif +!!$ +!!$ call psb_sp_reall(a,ia1_size,ia2_size,aspk_size,info) +!!$ if (info /= no_err) then +!!$ info=4010 +!!$ ch_err='psb_sp_reall' +!!$ call psb_errpush(info,name,a_err=ch_err) +!!$ goto 9999 +!!$ endif +!!$ +!!$ a%pl(:) = 0 +!!$ a%pr(:) = 0 if (debugwrt) then iout = 30+myrow