From 2c21b017b305ae6cb627f0ae58c79caeab5a77a3 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 18 Apr 2006 16:58:41 +0000 Subject: [PATCH] Fix complex version. --- src/prec/psb_zslu_impl.c | 28 ++- src/psblas/psb_zspsm.f90 | 2 +- src/serial/csr/dcsrsm.f | 2 +- src/serial/csr/zcsrsm.f | 90 +++++-- src/serial/csr/zsrmv.f | 516 +++++++++++++++++++-------------------- src/serial/csr/zsrsv.f | 242 +++++++++--------- src/serial/f77/zcssm.f | 2 + src/serial/psb_zcssv.f90 | 2 +- 8 files changed, 476 insertions(+), 408 deletions(-) diff --git a/src/prec/psb_zslu_impl.c b/src/prec/psb_zslu_impl.c index 0b383082..fa8398e4 100644 --- a/src/prec/psb_zslu_impl.c +++ b/src/prec/psb_zslu_impl.c @@ -85,7 +85,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #ifdef Have_SLU_ #include "zsp_defs.h" - +#undef Have_SLU_ #define HANDLE_SIZE 8 /* kind of integer to hold a pointer. Use int. This might need to be changed on 64-bit systems. */ @@ -102,7 +102,6 @@ typedef struct { int *perm_r; } factors_t; -#undef Have_SLU_ #else @@ -132,7 +131,7 @@ typedef struct { void psb_zslu_factor_(int *n, int *nnz, - double *values, int *rowind, int *colptr, + doublecomplex *values, int *rowind, int *colptr, #ifdef Have_SLU_ fptr *f_factors, /* a handle containing the address pointing to the factored matrices */ @@ -181,8 +180,8 @@ psb_zslu_factor_(int *n, int *nnz, for (i = 0; i < *nnz; ++i) --colptr[i]; for (i = 0; i <= *n; ++i) --rowind[i]; - dCreate_CompRow_Matrix(&A, *n, *n, *nnz, values, colptr, rowind, - SLU_NR, SLU_D, SLU_GE); + zCreate_CompRow_Matrix(&A, *n, *n, *nnz, values, colptr, rowind, + SLU_NC, SLU_Z, SLU_GE); L = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) ); U = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) ); if ( !(perm_r = intMalloc(*n)) ) ABORT("Malloc fails for perm_r[]."); @@ -205,13 +204,13 @@ psb_zslu_factor_(int *n, int *nnz, panel_size = sp_ienv(1); relax = sp_ienv(2); - dgstrf(&options, &AC, drop_tol, relax, panel_size, + zgstrf(&options, &AC, drop_tol, relax, panel_size, etree, NULL, 0, perm_c, perm_r, L, U, &stat, info); if ( *info == 0 ) { Lstore = (SCformat *) L->Store; Ustore = (NCformat *) U->Store; - dQuerySpace(L, U, &mem_usage); + zQuerySpace(L, U, &mem_usage); #if 0 printf("No of nonzeros in factor L = %d\n", Lstore->nnz); printf("No of nonzeros in factor U = %d\n", Ustore->nnz); @@ -223,7 +222,7 @@ psb_zslu_factor_(int *n, int *nnz, } else { printf("dgstrf() error returns INFO= %d\n", *info); if ( *info <= *n ) { /* factorization completes */ - dQuerySpace(L, U, &mem_usage); + zQuerySpace(L, U, &mem_usage); printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n", mem_usage.for_lu/1e6, mem_usage.total_needed/1e6, mem_usage.expansions); @@ -256,7 +255,7 @@ psb_zslu_factor_(int *n, int *nnz, void psb_zslu_solve_(int *itrans, int *n, int *nrhs, - double *b, int *ldb, + doublecomplex *b, int *ldb, #ifdef Have_SLU_ fptr *f_factors, /* a handle containing the address pointing to the factored matrices */ @@ -272,7 +271,7 @@ psb_zslu_solve_(int *itrans, int *n, int *nrhs, * */ #ifdef Have_SLU_ - SuperMatrix A, AC, B; + SuperMatrix B; SuperMatrix *L, *U; int *perm_r; /* row permutations from partial pivoting */ int *perm_c; /* column permutation vector */ @@ -306,9 +305,14 @@ psb_zslu_solve_(int *itrans, int *n, int *nrhs, perm_c = LUfactors->perm_c; perm_r = LUfactors->perm_r; - dCreate_Dense_Matrix(&B, *n, *nrhs, b, *ldb, SLU_DN, SLU_D, SLU_GE); + zCreate_Dense_Matrix(&B, *n, *nrhs, b, *ldb, SLU_DN, SLU_Z, SLU_GE); /* Solve the system A*X=B, overwriting B with X. */ - dgstrs (trans, L, U, perm_c, perm_r, &B, &stat, info); + zgstrs (trans, L, U, perm_c, perm_r, &B, &stat, info); + if (info != 0) { + if (B.Stype != SLU_DN) fprintf(stderr,"zgstrs error kind 1: SLU_DN\n"); + if (B.Dtype != SLU_Z) fprintf(stderr,"zgstrs error kind 2: SLU_Z\n"); + if (B.Mtype != SLU_GE) fprintf(stderr,"zgstrs error kind 3: SLU_GE\n"); + } Destroy_SuperMatrix_Store(&B); StatFree(&stat); diff --git a/src/psblas/psb_zspsm.f90 b/src/psblas/psb_zspsm.f90 index 28cf4a83..97abed25 100644 --- a/src/psblas/psb_zspsm.f90 +++ b/src/psblas/psb_zspsm.f90 @@ -262,7 +262,7 @@ subroutine psb_zspsm(alpha,a,x,beta,y,desc_a,info,& if(info.ne.0) then info = 4010 - ch_err='dcssm' + ch_err='zcssm' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if diff --git a/src/serial/csr/dcsrsm.f b/src/serial/csr/dcsrsm.f index 7fe4fc00..62251dc4 100644 --- a/src/serial/csr/dcsrsm.f +++ b/src/serial/csr/dcsrsm.f @@ -42,7 +42,7 @@ C PARAMETER (DEBUG=.FALSE.) C .. Local Arrays .. CHARACTER*20 NAME - INTEGER INT_VAL(5) + INTEGER INT_VAL(5), err_Act NAME = 'DCSRSM\0' IERROR = 0 diff --git a/src/serial/csr/zcsrsm.f b/src/serial/csr/zcsrsm.f index b8b05427..1aa22162 100644 --- a/src/serial/csr/zcsrsm.f +++ b/src/serial/csr/zcsrsm.f @@ -1,3 +1,33 @@ +C +C Parallel Sparse BLAS v2.0 +C (C) Copyright 2006 Salvatore Filippone University of Rome Tor Vergata +C Alfredo Buttari University of Rome Tor Vergata +C +C Redistribution and use in source and binary forms, with or without +C modification, are permitted provided that the following conditions +C are met: +C 1. Redistributions of source code must retain the above copyright +C notice, this list of conditions and the following disclaimer. +C 2. Redistributions in binary form must reproduce the above copyright +C notice, this list of conditions, and the following disclaimer in the +C documentation and/or other materials provided with the distribution. +C 3. The name of the PSBLAS group or the names of its contributors may +C not be used to endorse or promote products derived from this +C software without specific written permission. +C +C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +C BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +C CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +C POSSIBILITY OF SUCH DAMAGE. +C +C SUBROUTINE ZCSRSM(TRANST,M,N,UNITD,D,ALPHA,DESCRA,A,JA,IA, * B,LDB,BETA,C,LDC,WORK,LWORK) COMPLEX*16 ALPHA, BETA @@ -8,24 +38,31 @@ CHARACTER DESCRA*11 INTEGER I, K CHARACTER DIAG, UPLO - EXTERNAL XERBLA - IF ((ALPHA.NE.(1.D0, 0.D0)).OR. - + (BETA.NE.(0.D0, 0.D0))) THEN - CALL XERBLA('ZCSSM ',9) - RETURN - ENDIF + +C .. Local Arrays .. + CHARACTER*20 NAME + INTEGER INT_VAL(5),err_act + + NAME = 'DCSRSM\0' + IERROR = 0 + CALL FCPSB_ERRACTIONSAVE(ERR_ACT) + int_Val(1)=0 UPLO = '?' IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'U') UPLO = 'U' IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') UPLO = 'L' IF (UPLO.EQ.'?') THEN - CALL XERBLA('ZCSSM ',9) - RETURN + int_val(1) = 7 + IERROR=5 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 END IF IF (DESCRA(3:3).EQ.'N') DIAG = 'N' IF (DESCRA(3:3).EQ.'U') DIAG = 'U' IF(UNITD.EQ.'B') THEN - CALL XERBLA('ZCSSM ',9) - RETURN + IERROR=5 + int_val(1) = 4 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 ENDIF IF (UNITD.EQ.'R') THEN DO 40 I = 1, N @@ -35,9 +72,24 @@ 40 CONTINUE END IF - DO 60 I = 1, N - CALL ZSRSV(UPLO,TRANST,DIAG,M,A,JA,IA,B(1,I),C(1,I)) - 60 CONTINUE + if ((alpha.ne.(1.d0,0.d0)) .or.(beta .ne.(0.0d0,0.d0))) then + if (lwork .lt. m) then + int_val(1) = 17 + IERROR=5 + CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) + GOTO 9999 + END IF + DO I = 1, N + CALL ZSRSV(UPLO,TRANST,DIAG,M,A,JA,IA,B(1,I),work) + do k=1,m + c(k,i) = beta*c(k,i) + alpha*work(k) + enddo + enddo + else + DO I = 1, N + CALL ZSRSV(UPLO,TRANST,DIAG,M,A,JA,IA,B(1,I),c(1,i)) + enddo + endif IF (UNITD.EQ.'L') THEN DO 45 I = 1, N DO 25 K = 1, M @@ -45,7 +97,17 @@ 25 CONTINUE 45 CONTINUE END IF + + CALL FCPSB_ERRACTIONRESTORE(ERR_ACT) RETURN - END + 9999 CONTINUE + CALL FCPSB_ERRACTIONRESTORE(ERR_ACT) + + IF ( ERR_ACT .NE. 0 ) THEN + CALL FCPSB_SERROR() + RETURN + ENDIF + RETURN + END diff --git a/src/serial/csr/zsrmv.f b/src/serial/csr/zsrmv.f index 2188389c..b0e5e8c6 100644 --- a/src/serial/csr/zsrmv.f +++ b/src/serial/csr/zsrmv.f @@ -172,351 +172,351 @@ C .. Not simmetric matrix C .. Symmetric matrix upper or lower SYM = (TRANS.EQ.'L').OR.(TRANS.EQ.'U').OR. - + (TRANS.EQ.'M').OR.(TRANS.EQ.'V') + + (TRANS.EQ.'M').OR.(TRANS.EQ.'V') C IF (.NOT.(TRA.OR.COTRA)) THEN - NROWA = M - NCOLA = N + NROWA = M + NCOLA = N ELSE - NROWA = N - NCOLA = M - END IF !(....(CO)TRA) + NROWA = N + NCOLA = M + END IF !(....(CO)TRA) IF (ALPHA.EQ.ZERO) THEN - IF (BETA.EQ.ZERO) THEN - DO I = 1, M - Y(I) = ZERO - ENDDO - ELSE - DO 20 I = 1, M - Y(I) = BETA*Y(I) - 20 CONTINUE - ENDIF - RETURN + IF (BETA.EQ.ZERO) THEN + DO I = 1, M + Y(I) = ZERO + ENDDO + ELSE + DO 20 I = 1, M + Y(I) = BETA*Y(I) + 20 CONTINUE + ENDIF + RETURN END IF C IF (SYM) THEN - IF (UNI) THEN + IF (UNI) THEN C C ......Symmetric with unitary diagonal....... C ....OK!! C To be optimized - - IF (BETA.NE.ZERO) THEN - DO 40 I = 1, M + + IF (BETA.NE.ZERO) THEN + DO 40 I = 1, M C C Product for diagonal elements C - Y(I) = BETA*Y(I) + ALPHA*X(I) - 40 CONTINUE - ELSE - DO I = 1, M - Y(I) = ALPHA*X(I) - ENDDO - ENDIF + Y(I) = BETA*Y(I) + ALPHA*X(I) + 40 CONTINUE + ELSE + DO I = 1, M + Y(I) = ALPHA*X(I) + ENDDO + ENDIF C Product for other elements - IF ((TRANS.EQ.'L').OR.(TRANS.EQ.'U')) THEN - DO 80 I = 1, M - ACC = ZERO - DO 60 J = IA(I), IA(I+1) - 1 - K = JA(J) - Y(K) = Y(K) + ALPHA*AS(J)*X(I) - ACC = ACC + AS(J)*X(K) - 60 CONTINUE - Y(I) = Y(I) + ALPHA*ACC - 80 CONTINUE - ELSE ! Perform computations on conjug(A) - DO 82 I = 1, M - ACC = ZERO - DO 81 J = IA(I), IA(I+1) - 1 - K = JA(J) - Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) - ACC = ACC + CONJG(AS(J)) * X(K) - 81 CONTINUE - Y(I) = Y(I) + ALPHA*ACC - 82 CONTINUE - ENDIF -C - ELSE IF ( .NOT. UNI) THEN + IF ((TRANS.EQ.'L').OR.(TRANS.EQ.'U')) THEN + DO 80 I = 1, M + ACC = ZERO + DO 60 J = IA(I), IA(I+1) - 1 + K = JA(J) + Y(K) = Y(K) + ALPHA*AS(J)*X(I) + ACC = ACC + AS(J)*X(K) + 60 CONTINUE + Y(I) = Y(I) + ALPHA*ACC + 80 CONTINUE + ELSE ! Perform computations on conjug(A) + DO 82 I = 1, M + ACC = ZERO + DO 81 J = IA(I), IA(I+1) - 1 + K = JA(J) + Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) + ACC = ACC + CONJG(AS(J)) * X(K) + 81 CONTINUE + Y(I) = Y(I) + ALPHA*ACC + 82 CONTINUE + ENDIF +C + ELSE IF ( .NOT. UNI) THEN C C Check if matrix is lower or upper C - IF ((TRANS.EQ.'L').OR.(TRANS.EQ.'M')) THEN + IF ((TRANS.EQ.'L').OR.(TRANS.EQ.'M')) THEN C C LOWER CASE: diagonal element is the last element of row C ....OK! - IF (BETA.NE.ZERO) THEN - DO 100 I = 1, M - Y(I) = BETA*Y(I) - 100 CONTINUE - ELSE - DO I = 1, M - Y(I) = ZERO - ENDDO - ENDIF + IF (BETA.NE.ZERO) THEN + DO 100 I = 1, M + Y(I) = BETA*Y(I) + 100 CONTINUE + ELSE + DO I = 1, M + Y(I) = ZERO + ENDDO + ENDIF - IF (TRANS.EQ.'L') THEN - DO 140 I = 1, M - ACC = ZERO - DO 120 J = IA(I), IA(I+1) - 1 ! it was -2 - K = JA(J) + IF (TRANS.EQ.'L') THEN + DO 140 I = 1, M + ACC = ZERO + DO 120 J = IA(I), IA(I+1) - 1 ! it was -2 + K = JA(J) C C To be optimized C - IF (K.NE.I) THEN !for symmetric element - Y(K) = Y(K) + ALPHA*AS(J)*X(I) - ENDIF - ACC = ACC + AS(J)*X(K) - 120 CONTINUE + IF (K.NE.I) THEN !for symmetric element + Y(K) = Y(K) + ALPHA*AS(J)*X(I) + ENDIF + ACC = ACC + AS(J)*X(K) + 120 CONTINUE - Y(I) = Y(I) + ALPHA*ACC - 140 CONTINUE - ELSE ! Perform computations on conjug(A) - DO 142 I = 1, M - ACC = ZERO - DO 141 J = IA(I), IA(I+1) - 1 ! it was -2 - K = JA(J) + Y(I) = Y(I) + ALPHA*ACC + 140 CONTINUE + ELSE ! Perform computations on conjug(A) + DO 142 I = 1, M + ACC = ZERO + DO 141 J = IA(I), IA(I+1) - 1 ! it was -2 + K = JA(J) C C To be optimized C - IF (K.NE.I) THEN !for symmetric element - Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) - ENDIF - ACC = ACC + CONJG(AS(J)) * X(K) - 141 CONTINUE - Y(I) = Y(I) + ALPHA * ACC - 142 CONTINUE + IF (K.NE.I) THEN !for symmetric element + Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) + ENDIF + ACC = ACC + CONJG(AS(J)) * X(K) + 141 CONTINUE + Y(I) = Y(I) + ALPHA * ACC + 142 CONTINUE - ENDIF + ENDIF - ELSE ! ....Trans<>L, M + ELSE ! ....Trans<>L, M C C UPPER CASE C ....OK!! C - IF (BETA.NE.ZERO) THEN - DO 160 I = 1, M - Y(I) = BETA*Y(I) - 160 CONTINUE - ELSE - DO I = 1, M - Y(I) = ZERO - ENDDO - ENDIF - IF (TRANS.EQ.'U') THEN - DO 200 I = 1, M - ACC = ZERO - DO 180 J = IA(I) , IA(I+1) - 1 ! removed +1 - K = JA(J) + IF (BETA.NE.ZERO) THEN + DO 160 I = 1, M + Y(I) = BETA*Y(I) + 160 CONTINUE + ELSE + DO I = 1, M + Y(I) = ZERO + ENDDO + ENDIF + IF (TRANS.EQ.'U') THEN + DO 200 I = 1, M + ACC = ZERO + DO 180 J = IA(I) , IA(I+1) - 1 ! removed +1 + K = JA(J) C C To be optimized C - IF(K.NE.I) THEN - Y(K) = Y(K) + ALPHA*AS(J)*X(I) - ENDIF - ACC = ACC + AS(J)*X(K) - 180 CONTINUE - Y(I) = Y(I) + ALPHA*ACC - 200 CONTINUE - ELSE ! Perform computations on conjug(A) - DO 202 I = 1, M - ACC = ZERO - DO 201 J = IA(I) , IA(I+1) - 1 ! removed +1 - K = JA(J) + IF(K.NE.I) THEN + Y(K) = Y(K) + ALPHA*AS(J)*X(I) + ENDIF + ACC = ACC + AS(J)*X(K) + 180 CONTINUE + Y(I) = Y(I) + ALPHA*ACC + 200 CONTINUE + ELSE ! Perform computations on conjug(A) + DO 202 I = 1, M + ACC = ZERO + DO 201 J = IA(I) , IA(I+1) - 1 ! removed +1 + K = JA(J) C C To be optimized C - IF(K.NE.I) THEN - Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) - ENDIF - ACC = ACC + CONJG(AS(J)) * X(K) - 201 CONTINUE - Y(I) = Y(I) + ALPHA*ACC - 202 CONTINUE - ENDIF - END IF ! ......TRANS=='L' - END IF ! ......Not UNI -C - ELSE IF ((.NOT.TRA).AND.(.NOT.COTRA)) THEN !......NOT SYM + IF(K.NE.I) THEN + Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) + ENDIF + ACC = ACC + CONJG(AS(J)) * X(K) + 201 CONTINUE + Y(I) = Y(I) + ALPHA*ACC + 202 CONTINUE + ENDIF + END IF ! ......TRANS=='L' + END IF ! ......Not UNI +C + ELSE IF ((.NOT.TRA).AND.(.NOT.COTRA)) THEN !......NOT SYM - IF ( .NOT. UNI) THEN + IF ( .NOT. UNI) THEN C C .......General Not Unit, No Traspose C - IF (BETA.NE.ZERO) THEN - DO 240 I = 1, M - ACC = ZERO - DO 220 J = IA(I), IA(I+1) - 1 - ACC = ACC + AS(J)*X(JA(J)) - 220 CONTINUE - Y(I) = ALPHA*ACC + BETA*Y(I) - 240 CONTINUE - ELSE - DO I = 1, M - ACC = ZERO - DO J = IA(I), IA(I+1) - 1 - ACC = ACC + AS(J)*X(JA(J)) - ENDDO - Y(I) = ALPHA*ACC - ENDDO - ENDIF -C - ELSE IF (UNI) THEN -C - IF (BETA.NE.ZERO) THEN - DO 280 I = 1, M - ACC = ZERO - DO 260 J = IA(I), IA(I+1) - 1 - ACC = ACC + AS(J)*X(JA(J)) - 260 CONTINUE - Y(I) = ALPHA*(ACC+X(I)) + BETA*Y(I) - 280 CONTINUE - ELSE !(BETA.EQ.ZERO) - DO I = 1, M - ACC = ZERO - DO J = IA(I), IA(I+1) - 1 - ACC = ACC + AS(J)*X(JA(J)) - ENDDO - Y(I) = ALPHA*(ACC+X(I)) - ENDDO - ENDIF - END IF !....End Testing on UNI + IF (BETA.NE.ZERO) THEN + DO 240 I = 1, M + ACC = ZERO + DO 220 J = IA(I), IA(I+1) - 1 + ACC = ACC + AS(J)*X(JA(J)) + 220 CONTINUE + Y(I) = ALPHA*ACC + BETA*Y(I) + 240 CONTINUE + ELSE + DO I = 1, M + ACC = ZERO + DO J = IA(I), IA(I+1) - 1 + ACC = ACC + AS(J)*X(JA(J)) + ENDDO + Y(I) = ALPHA*ACC + ENDDO + ENDIF +C + ELSE IF (UNI) THEN +C + IF (BETA.NE.ZERO) THEN + DO 280 I = 1, M + ACC = ZERO + DO 260 J = IA(I), IA(I+1) - 1 + ACC = ACC + AS(J)*X(JA(J)) + 260 CONTINUE + Y(I) = ALPHA*(ACC+X(I)) + BETA*Y(I) + 280 CONTINUE + ELSE !(BETA.EQ.ZERO) + DO I = 1, M + ACC = ZERO + DO J = IA(I), IA(I+1) - 1 + ACC = ACC + AS(J)*X(JA(J)) + ENDDO + Y(I) = ALPHA*(ACC+X(I)) + ENDDO + ENDIF + END IF !....End Testing on UNI C - ELSE IF (TRA) THEN !....Else on SYM (swapped M and N) + ELSE IF (TRA) THEN !....Else on SYM (swapped M and N) C - IF ( .NOT. UNI) THEN + IF ( .NOT. UNI) THEN C - IF (BETA.NE.ZERO) THEN - DO 300 I = 1, M - Y(I) = BETA*Y(I) - 300 CONTINUE - ELSE !(BETA.EQ.ZERO) - DO I = 1, M - Y(I) = ZERO - ENDDO - ENDIF + IF (BETA.NE.ZERO) THEN + DO 300 I = 1, M + Y(I) = BETA*Y(I) + 300 CONTINUE + ELSE !(BETA.EQ.ZERO) + DO I = 1, M + Y(I) = ZERO + ENDDO + ENDIF C - ELSE IF (UNI) THEN + ELSE IF (UNI) THEN C - IF (BETA.NE.ZERO) THEN - DO 320 I = 1, M - Y(I) = BETA*Y(I) + ALPHA*X(I) - 320 CONTINUE - ELSE !(BETA.EQ.ZERO) - DO I = 1, M - Y(I) = ALPHA*X(I) - ENDDO - ENDIF + IF (BETA.NE.ZERO) THEN + DO 320 I = 1, M + Y(I) = BETA*Y(I) + ALPHA*X(I) + 320 CONTINUE + ELSE !(BETA.EQ.ZERO) + DO I = 1, M + Y(I) = ALPHA*X(I) + ENDDO + ENDIF C - END IF !....UNI + END IF !....UNI C - IF (ALPHA.EQ.ONE) THEN + IF (ALPHA.EQ.ONE) THEN C - DO 360 I = 1, N - DO 340 J = IA(I), IA(I+1) - 1 - K = JA(J) - Y(K) = Y(K) + AS(J)*X(I) - 340 CONTINUE - 360 CONTINUE + DO 360 I = 1, N + DO 340 J = IA(I), IA(I+1) - 1 + K = JA(J) + Y(K) = Y(K) + AS(J)*X(I) + 340 CONTINUE + 360 CONTINUE C - ELSE IF (ALPHA.EQ.-ONE) THEN + ELSE IF (ALPHA.EQ.-ONE) THEN C - DO 400 I = 1, n - DO 380 J = IA(I), IA(I+1) - 1 - K = JA(J) - Y(K) = Y(K) - AS(J)*X(I) - 380 CONTINUE - 400 CONTINUE + DO 400 I = 1, n + DO 380 J = IA(I), IA(I+1) - 1 + K = JA(J) + Y(K) = Y(K) - AS(J)*X(I) + 380 CONTINUE + 400 CONTINUE C - ELSE !.....Else on TRA + ELSE !.....Else on TRA C C This work array is used for efficiency C - DO 420 I = 1, N - WORK(I) = ALPHA*X(I) - 420 CONTINUE + DO 420 I = 1, N + WORK(I) = ALPHA*X(I) + 420 CONTINUE C - DO 460 I = 1, n - DO 440 J = IA(I), IA(I+1) - 1 - K = JA(J) - Y(K) = Y(K) + AS(J)*WORK(I) - 440 CONTINUE - 460 CONTINUE + DO 460 I = 1, n + DO 440 J = IA(I), IA(I+1) - 1 + K = JA(J) + Y(K) = Y(K) + AS(J)*WORK(I) + 440 CONTINUE + 460 CONTINUE C - END IF !.....End testing on ALPHA + END IF !.....End testing on ALPHA - ELSE IF (COTRA) THEN !....Else on SYM (swapped M and N) + ELSE IF (COTRA) THEN !....Else on SYM (swapped M and N) C - IF ( .NOT. UNI) THEN + IF ( .NOT. UNI) THEN C - IF (BETA.NE.ZERO) THEN - DO 500 I = 1, M - Y(I) = BETA*Y(I) - 500 CONTINUE - ELSE !(BETA.EQ.ZERO) - DO I = 1, M - Y(I) = ZERO - ENDDO - ENDIF + IF (BETA.NE.ZERO) THEN + DO 500 I = 1, M + Y(I) = BETA*Y(I) + 500 CONTINUE + ELSE !(BETA.EQ.ZERO) + DO I = 1, M + Y(I) = ZERO + ENDDO + ENDIF C - ELSE IF (UNI) THEN + ELSE IF (UNI) THEN C - IF (BETA.NE.ZERO) THEN - DO 520 I = 1, M - Y(I) = BETA*Y(I) + ALPHA*X(I) - 520 CONTINUE - ELSE !(BETA.EQ.ZERO) - DO I = 1, M - Y(I) = ALPHA*X(I) - ENDDO - ENDIF + IF (BETA.NE.ZERO) THEN + DO 520 I = 1, M + Y(I) = BETA*Y(I) + ALPHA*X(I) + 520 CONTINUE + ELSE !(BETA.EQ.ZERO) + DO I = 1, M + Y(I) = ALPHA*X(I) + ENDDO + ENDIF C - END IF !....UNI + END IF !....UNI C - IF (ALPHA.EQ.ONE) THEN + IF (ALPHA.EQ.ONE) THEN C - DO 560 I = 1, N - DO 540 J = IA(I), IA(I+1) - 1 - K = JA(J) - Y(K) = Y(K) + CONJG(AS(J))*X(I) - 540 CONTINUE - 560 CONTINUE + DO 560 I = 1, N + DO 540 J = IA(I), IA(I+1) - 1 + K = JA(J) + Y(K) = Y(K) + CONJG(AS(J))*X(I) + 540 CONTINUE + 560 CONTINUE C - ELSE IF (ALPHA.EQ.-ONE) THEN + ELSE IF (ALPHA.EQ.-ONE) THEN C - DO 600 I = 1, n - DO 580 J = IA(I), IA(I+1) - 1 - K = JA(J) - Y(K) = Y(K) - CONJG(AS(J))*X(I) - 580 CONTINUE - 600 CONTINUE + DO 600 I = 1, n + DO 580 J = IA(I), IA(I+1) - 1 + K = JA(J) + Y(K) = Y(K) - CONJG(AS(J))*X(I) + 580 CONTINUE + 600 CONTINUE C - ELSE !.....Else on TRA + ELSE !.....Else on TRA C C This work array is used for efficiency C - DO 620 I = 1, N - WORK(I) = ALPHA*X(I) - 620 CONTINUE + DO 620 I = 1, N + WORK(I) = ALPHA*X(I) + 620 CONTINUE C - DO 660 I = 1, n - DO 640 J = IA(I), IA(I+1) - 1 - K = JA(J) - Y(K) = Y(K) + CONJG(AS(J))*WORK(I) - 640 CONTINUE - 660 CONTINUE + DO 660 I = 1, n + DO 640 J = IA(I), IA(I+1) - 1 + K = JA(J) + Y(K) = Y(K) + CONJG(AS(J))*WORK(I) + 640 CONTINUE + 660 CONTINUE C - END IF !.....End testing on ALPHA + END IF !.....End testing on ALPHA - END IF !.....End testing on SYM + END IF !.....End testing on SYM C RETURN C diff --git a/src/serial/csr/zsrsv.f b/src/serial/csr/zsrsv.f index 32b484bc..cf0e5158 100644 --- a/src/serial/csr/zsrsv.f +++ b/src/serial/csr/zsrsv.f @@ -13,129 +13,129 @@ COTRA = (TRANS.EQ.'C') LOW = (UPLO.EQ.'L') IF ((.NOT.TRA).AND.(.NOT.COTRA)) THEN - IF (LOW) THEN - IF (.NOT.UNI) THEN - DO 40 I = 1, N - ACC = ZERO - DO 20 J = IA(I), IA(I+1) - 2 - ACC = ACC + AS(J)*X(JA(J)) - 20 CONTINUE - X(I) = (B(I)-ACC)/AS(IA(I+1)-1) - 40 CONTINUE - ELSE IF (UNI) THEN - DO 80 I = 1, N - ACC = ZERO - DO 60 J = IA(I), IA(I+1) - 1 - ACC = ACC + AS(J)*X(JA(J)) - 60 CONTINUE - X(I) = B(I) - ACC - 80 CONTINUE - END IF - ELSE IF (.NOT.LOW) THEN - IF (.NOT.UNI) THEN - DO 120 I = N, 1, -1 - ACC = ZERO - DO 100 J = IA(I) + 1, IA(I+1) - 1 - ACC = ACC + AS(J)*X(JA(J)) - 100 CONTINUE - X(I) = (B(I)-ACC)/AS(IA(I)) - 120 CONTINUE - ELSE IF (UNI) THEN - DO 160 I = N, 1, -1 - ACC = ZERO - DO 140 J = IA(I), IA(I+1) - 1 - ACC = ACC + AS(J)*X(JA(J)) - 140 CONTINUE - X(I) = B(I) - ACC - 160 CONTINUE - END IF - END IF + IF (LOW) THEN + IF (.NOT.UNI) THEN + DO 40 I = 1, N + ACC = ZERO + DO 20 J = IA(I), IA(I+1) - 2 + ACC = ACC + AS(J)*X(JA(J)) + 20 CONTINUE + X(I) = (B(I)-ACC)/AS(IA(I+1)-1) + 40 CONTINUE + ELSE IF (UNI) THEN + DO 80 I = 1, N + ACC = ZERO + DO 60 J = IA(I), IA(I+1) - 1 + ACC = ACC + AS(J)*X(JA(J)) + 60 CONTINUE + X(I) = B(I) - ACC + 80 CONTINUE + END IF + ELSE IF (.NOT.LOW) THEN + IF (.NOT.UNI) THEN + DO 120 I = N, 1, -1 + ACC = ZERO + DO 100 J = IA(I) + 1, IA(I+1) - 1 + ACC = ACC + AS(J)*X(JA(J)) + 100 CONTINUE + X(I) = (B(I)-ACC)/AS(IA(I)) + 120 CONTINUE + ELSE IF (UNI) THEN + DO 160 I = N, 1, -1 + ACC = ZERO + DO 140 J = IA(I), IA(I+1) - 1 + ACC = ACC + AS(J)*X(JA(J)) + 140 CONTINUE + X(I) = B(I) - ACC + 160 CONTINUE + END IF + END IF ELSE IF (TRA) THEN - DO 180 I = 1, N - X(I) = B(I) - 180 CONTINUE - IF (LOW) THEN - IF (.NOT.UNI) THEN - DO 220 I = N, 1, -1 - X(I) = X(I)/AS(IA(I+1)-1) - ACC = X(I) - DO 200 J = IA(I), IA(I+1) - 2 - K = JA(J) - X(K) = X(K) - AS(J)*ACC - 200 CONTINUE - 220 CONTINUE - ELSE IF (UNI) THEN - DO 260 I = N, 1, -1 - ACC = X(I) - DO 240 J = IA(I), IA(I+1) - 1 - K = JA(J) - X(K) = X(K) - AS(J)*ACC - 240 CONTINUE - 260 CONTINUE - END IF - ELSE IF (.NOT.LOW) THEN - IF (.NOT.UNI) THEN - DO 300 I = 1, N - X(I) = X(I)/AS(IA(I)) - ACC = X(I) - DO 280 J = IA(I) + 1, IA(I+1) - 1 - K = JA(J) - X(K) = X(K) - AS(J)*ACC - 280 CONTINUE - 300 CONTINUE - ELSE IF (UNI) THEN - DO 340 I = 1, N - ACC = X(I) - DO 320 J = IA(I), IA(I+1) - 1 - K = JA(J) - X(K) = X(K) - AS(J)*ACC - 320 CONTINUE - 340 CONTINUE - END IF - END IF + DO 180 I = 1, N + X(I) = B(I) + 180 CONTINUE + IF (LOW) THEN + IF (.NOT.UNI) THEN + DO 220 I = N, 1, -1 + X(I) = X(I)/AS(IA(I+1)-1) + ACC = X(I) + DO 200 J = IA(I), IA(I+1) - 2 + K = JA(J) + X(K) = X(K) - AS(J)*ACC + 200 CONTINUE + 220 CONTINUE + ELSE IF (UNI) THEN + DO 260 I = N, 1, -1 + ACC = X(I) + DO 240 J = IA(I), IA(I+1) - 1 + K = JA(J) + X(K) = X(K) - AS(J)*ACC + 240 CONTINUE + 260 CONTINUE + END IF + ELSE IF (.NOT.LOW) THEN + IF (.NOT.UNI) THEN + DO 300 I = 1, N + X(I) = X(I)/AS(IA(I)) + ACC = X(I) + DO 280 J = IA(I) + 1, IA(I+1) - 1 + K = JA(J) + X(K) = X(K) - AS(J)*ACC + 280 CONTINUE + 300 CONTINUE + ELSE IF (UNI) THEN + DO 340 I = 1, N + ACC = X(I) + DO 320 J = IA(I), IA(I+1) - 1 + K = JA(J) + X(K) = X(K) - AS(J)*ACC + 320 CONTINUE + 340 CONTINUE + END IF + END IF ELSE IF (COTRA) THEN - DO 580 I = 1, N - X(I) = B(I) - 580 CONTINUE - IF (LOW) THEN - IF (.NOT.UNI) THEN - DO 620 I = N, 1, -1 - X(I) = X(I)/CONJG(AS(IA(I+1)-1)) - ACC = X(I) - DO 600 J = IA(I), IA(I+1) - 2 - K = JA(J) - X(K) = X(K) - CONJG(AS(J))*ACC - 600 CONTINUE - 620 CONTINUE - ELSE IF (UNI) THEN - DO 660 I = N, 1, -1 - ACC = X(I) - DO 640 J = IA(I), IA(I+1) - 1 - K = JA(J) - X(K) = X(K) - CONJG(AS(J))*ACC - 640 CONTINUE - 660 CONTINUE - END IF - ELSE IF (.NOT.LOW) THEN - IF (.NOT.UNI) THEN - DO 700 I = 1, N - X(I) = X(I)/CONJG(AS(IA(I))) - ACC = X(I) - DO 680 J = IA(I) + 1, IA(I+1) - 1 - K = JA(J) - X(K) = X(K) - CONJG(AS(J))*ACC - 680 CONTINUE - 700 CONTINUE - ELSE IF (UNI) THEN - DO 740 I = 1, N - ACC = X(I) - DO 720 J = IA(I), IA(I+1) - 1 - K = JA(J) - X(K) = X(K) - CONJG(AS(J))*ACC - 720 CONTINUE - 740 CONTINUE - END IF - END IF + DO 580 I = 1, N + X(I) = B(I) + 580 CONTINUE + IF (LOW) THEN + IF (.NOT.UNI) THEN + DO 620 I = N, 1, -1 + X(I) = X(I)/CONJG(AS(IA(I+1)-1)) + ACC = X(I) + DO 600 J = IA(I), IA(I+1) - 2 + K = JA(J) + X(K) = X(K) - CONJG(AS(J))*ACC + 600 CONTINUE + 620 CONTINUE + ELSE IF (UNI) THEN + DO 660 I = N, 1, -1 + ACC = X(I) + DO 640 J = IA(I), IA(I+1) - 1 + K = JA(J) + X(K) = X(K) - CONJG(AS(J))*ACC + 640 CONTINUE + 660 CONTINUE + END IF + ELSE IF (.NOT.LOW) THEN + IF (.NOT.UNI) THEN + DO 700 I = 1, N + X(I) = X(I)/CONJG(AS(IA(I))) + ACC = X(I) + DO 680 J = IA(I) + 1, IA(I+1) - 1 + K = JA(J) + X(K) = X(K) - CONJG(AS(J))*ACC + 680 CONTINUE + 700 CONTINUE + ELSE IF (UNI) THEN + DO 740 I = 1, N + ACC = X(I) + DO 720 J = IA(I), IA(I+1) - 1 + K = JA(J) + X(K) = X(K) - CONJG(AS(J))*ACC + 720 CONTINUE + 740 CONTINUE + END IF + END IF END IF RETURN END diff --git a/src/serial/f77/zcssm.f b/src/serial/f77/zcssm.f index 00078940..a7d33ce5 100644 --- a/src/serial/f77/zcssm.f +++ b/src/serial/f77/zcssm.f @@ -231,6 +231,7 @@ C C Error handling C IF(IERROR.NE.0) THEN + write(0,*) 'zcssm ierror',ierror CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) GOTO 9999 @@ -335,6 +336,7 @@ C RETURN 9999 CONTINUE CALL FCPSB_ERRACTIONRESTORE(ERR_ACT) + IF ( ERR_ACT .NE. 0 ) THEN CALL FCPSB_SERROR() diff --git a/src/serial/psb_zcssv.f90 b/src/serial/psb_zcssv.f90 index 1c85bfbe..6bf16654 100644 --- a/src/serial/psb_zcssv.f90 +++ b/src/serial/psb_zcssv.f90 @@ -80,7 +80,7 @@ subroutine psb_zcssv(alpha,t,b,beta,c,info,trans,unitd,d) call zcssm(lt,m,n,alpha,lu,ddl,& & t%pl,t%fida,t%descra,t%aspk,t%ia1,t%ia2,t%infoa,t%pr,& & b,lb,beta,c,lc,work,iwsz,info) - + if (.not.present(d)) then deallocate(ddl) endif