fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fix complex version.

Browse Source
psblas3-type-indexed
Salvatore Filippone 19 years ago
parent 8ca0d6fa93
commit 2c21b017b3
8 changed files with 476 additions and 408 deletions
Show Stats Download Patch File Download Diff File

28
src/prec/psb_zslu_impl.c
Unescape Escape View File

@ -85,7 +85,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#ifdef Have_SLU_ #ifdef Have_SLU_
#include "zsp_defs.h" #include "zsp_defs.h"
#undef Have_SLU_
#define HANDLE_SIZE 8 #define HANDLE_SIZE 8
/* kind of integer to hold a pointer. Use int. /* kind of integer to hold a pointer. Use int.
This might need to be changed on 64-bit systems. */ This might need to be changed on 64-bit systems. */
@ -102,7 +102,6 @@ typedef struct {
int *perm_r; int *perm_r;
} factors_t; } factors_t;
#undef Have_SLU_
#else #else
@ -132,7 +131,7 @@ typedef struct {
void void
psb_zslu_factor_(int *n, int *nnz, psb_zslu_factor_(int *n, int *nnz,
double *values, int *rowind, int *colptr, doublecomplex *values, int *rowind, int *colptr,
#ifdef Have_SLU_ #ifdef Have_SLU_
fptr *f_factors, /* a handle containing the address fptr *f_factors, /* a handle containing the address
pointing to the factored matrices */ 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 < *nnz; ++i) --colptr[i];
for (i = 0; i <= *n; ++i) --rowind[i]; for (i = 0; i <= *n; ++i) --rowind[i];
dCreate_CompRow_Matrix(&A, *n, *n, *nnz, values, colptr, rowind, zCreate_CompRow_Matrix(&A, *n, *n, *nnz, values, colptr, rowind,
SLU_NR, SLU_D, SLU_GE); SLU_NC, SLU_Z, SLU_GE);
L = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) ); L = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
U = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) ); U = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
if ( !(perm_r = intMalloc(*n)) ) ABORT("Malloc fails for perm_r[]."); 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); panel_size = sp_ienv(1);
relax = sp_ienv(2); 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); etree, NULL, 0, perm_c, perm_r, L, U, &stat, info);
if ( *info == 0 ) { if ( *info == 0 ) {
Lstore = (SCformat *) L->Store; Lstore = (SCformat *) L->Store;
Ustore = (NCformat *) U->Store; Ustore = (NCformat *) U->Store;
dQuerySpace(L, U, &mem_usage); zQuerySpace(L, U, &mem_usage);
#if 0 #if 0
printf("No of nonzeros in factor L = %d\n", Lstore->nnz); printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
printf("No of nonzeros in factor U = %d\n", Ustore->nnz); printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
@ -223,7 +222,7 @@ psb_zslu_factor_(int *n, int *nnz,
} else { } else {
printf("dgstrf() error returns INFO= %d\n", *info); printf("dgstrf() error returns INFO= %d\n", *info);
if ( *info <= *n ) { /* factorization completes */ 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", printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6, mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions); mem_usage.expansions);
@ -256,7 +255,7 @@ psb_zslu_factor_(int *n, int *nnz,
void void
psb_zslu_solve_(int *itrans, int *n, int *nrhs, psb_zslu_solve_(int *itrans, int *n, int *nrhs,
double *b, int *ldb, doublecomplex *b, int *ldb,
#ifdef Have_SLU_ #ifdef Have_SLU_
fptr *f_factors, /* a handle containing the address fptr *f_factors, /* a handle containing the address
pointing to the factored matrices */ pointing to the factored matrices */
@ -272,7 +271,7 @@ psb_zslu_solve_(int *itrans, int *n, int *nrhs,
* *
*/ */
#ifdef Have_SLU_ #ifdef Have_SLU_
SuperMatrix A, AC, B; SuperMatrix B;
SuperMatrix *L, *U; SuperMatrix *L, *U;
int *perm_r; /* row permutations from partial pivoting */ int *perm_r; /* row permutations from partial pivoting */
int *perm_c; /* column permutation vector */ 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_c = LUfactors->perm_c;
perm_r = LUfactors->perm_r; 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. */ /* 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); Destroy_SuperMatrix_Store(&B);
StatFree(&stat); StatFree(&stat);

2
src/psblas/psb_zspsm.f90
Unescape Escape View File

@ -262,7 +262,7 @@ subroutine psb_zspsm(alpha,a,x,beta,y,desc_a,info,&
if(info.ne.0) then if(info.ne.0) then
info = 4010 info = 4010
ch_err='dcssm' ch_err='zcssm'
call psb_errpush(info,name,a_err=ch_err) call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if end if

2
src/serial/csr/dcsrsm.f
Unescape Escape View File

@ -42,7 +42,7 @@ C
PARAMETER (DEBUG=.FALSE.) PARAMETER (DEBUG=.FALSE.)
C .. Local Arrays .. C .. Local Arrays ..
CHARACTER*20 NAME CHARACTER*20 NAME
INTEGER INT_VAL(5) INTEGER INT_VAL(5), err_Act
NAME = 'DCSRSM\0' NAME = 'DCSRSM\0'
IERROR = 0 IERROR = 0

90
src/serial/csr/zcsrsm.f
Unescape Escape View File

@ -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, SUBROUTINE ZCSRSM(TRANST,M,N,UNITD,D,ALPHA,DESCRA,A,JA,IA,
* B,LDB,BETA,C,LDC,WORK,LWORK) * B,LDB,BETA,C,LDC,WORK,LWORK)
COMPLEX*16 ALPHA, BETA COMPLEX*16 ALPHA, BETA
@ -8,24 +38,31 @@
CHARACTER DESCRA*11 CHARACTER DESCRA*11
INTEGER I, K INTEGER I, K
CHARACTER DIAG, UPLO CHARACTER DIAG, UPLO
EXTERNAL XERBLA
IF ((ALPHA.NE.(1.D0, 0.D0)).OR. C .. Local Arrays ..
+ (BETA.NE.(0.D0, 0.D0))) THEN CHARACTER*20 NAME
CALL XERBLA('ZCSSM ',9) INTEGER INT_VAL(5),err_act
RETURN
ENDIF NAME = 'DCSRSM\0'
IERROR = 0
CALL FCPSB_ERRACTIONSAVE(ERR_ACT)
int_Val(1)=0
UPLO = '?' 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.'U') UPLO = 'U'
IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') UPLO = 'L' IF (DESCRA(1:1).EQ.'T' .AND. DESCRA(2:2).EQ.'L') UPLO = 'L'
IF (UPLO.EQ.'?') THEN IF (UPLO.EQ.'?') THEN
CALL XERBLA('ZCSSM ',9) int_val(1) = 7
RETURN IERROR=5
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
END IF END IF
IF (DESCRA(3:3).EQ.'N') DIAG = 'N' IF (DESCRA(3:3).EQ.'N') DIAG = 'N'
IF (DESCRA(3:3).EQ.'U') DIAG = 'U' IF (DESCRA(3:3).EQ.'U') DIAG = 'U'
IF(UNITD.EQ.'B') THEN IF(UNITD.EQ.'B') THEN
CALL XERBLA('ZCSSM ',9) IERROR=5
RETURN int_val(1) = 4
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999
ENDIF ENDIF
IF (UNITD.EQ.'R') THEN IF (UNITD.EQ.'R') THEN
DO 40 I = 1, N DO 40 I = 1, N
@ -35,9 +72,24 @@
40 CONTINUE 40 CONTINUE
END IF END IF
DO 60 I = 1, N if ((alpha.ne.(1.d0,0.d0)) .or.(beta .ne.(0.0d0,0.d0))) then
CALL ZSRSV(UPLO,TRANST,DIAG,M,A,JA,IA,B(1,I),C(1,I)) if (lwork .lt. m) then
60 CONTINUE 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 IF (UNITD.EQ.'L') THEN
DO 45 I = 1, N DO 45 I = 1, N
DO 25 K = 1, M DO 25 K = 1, M
@ -45,7 +97,17 @@
25 CONTINUE 25 CONTINUE
45 CONTINUE 45 CONTINUE
END IF END IF
CALL FCPSB_ERRACTIONRESTORE(ERR_ACT)
RETURN RETURN
END
9999 CONTINUE
CALL FCPSB_ERRACTIONRESTORE(ERR_ACT)
IF ( ERR_ACT .NE. 0 ) THEN
CALL FCPSB_SERROR()
RETURN
ENDIF
RETURN
END

516
src/serial/csr/zsrmv.f
Unescape Escape View File

@ -172,351 +172,351 @@ C .. Not simmetric matrix
C .. Symmetric matrix upper or lower C .. Symmetric matrix upper or lower
SYM = (TRANS.EQ.'L').OR.(TRANS.EQ.'U').OR. SYM = (TRANS.EQ.'L').OR.(TRANS.EQ.'U').OR.
+ (TRANS.EQ.'M').OR.(TRANS.EQ.'V') + (TRANS.EQ.'M').OR.(TRANS.EQ.'V')
C C
IF (.NOT.(TRA.OR.COTRA)) THEN IF (.NOT.(TRA.OR.COTRA)) THEN
NROWA = M NROWA = M
NCOLA = N NCOLA = N
ELSE ELSE
NROWA = N NROWA = N
NCOLA = M NCOLA = M
END IF !(....(CO)TRA) END IF !(....(CO)TRA)
IF (ALPHA.EQ.ZERO) THEN IF (ALPHA.EQ.ZERO) THEN
IF (BETA.EQ.ZERO) THEN IF (BETA.EQ.ZERO) THEN
DO I = 1, M DO I = 1, M
Y(I) = ZERO Y(I) = ZERO
ENDDO ENDDO
ELSE ELSE
DO 20 I = 1, M DO 20 I = 1, M
Y(I) = BETA*Y(I) Y(I) = BETA*Y(I)
20 CONTINUE 20 CONTINUE
ENDIF ENDIF
RETURN RETURN
END IF END IF
C C
IF (SYM) THEN IF (SYM) THEN
IF (UNI) THEN IF (UNI) THEN
C C
C ......Symmetric with unitary diagonal....... C ......Symmetric with unitary diagonal.......
C ....OK!! C ....OK!!
C To be optimized C To be optimized
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 40 I = 1, M DO 40 I = 1, M
C C
C Product for diagonal elements C Product for diagonal elements
C C
Y(I) = BETA*Y(I) + ALPHA*X(I) Y(I) = BETA*Y(I) + ALPHA*X(I)
40 CONTINUE 40 CONTINUE
ELSE ELSE
DO I = 1, M DO I = 1, M
Y(I) = ALPHA*X(I) Y(I) = ALPHA*X(I)
ENDDO ENDDO
ENDIF ENDIF
C Product for other elements C Product for other elements
IF ((TRANS.EQ.'L').OR.(TRANS.EQ.'U')) THEN IF ((TRANS.EQ.'L').OR.(TRANS.EQ.'U')) THEN
DO 80 I = 1, M DO 80 I = 1, M
ACC = ZERO ACC = ZERO
DO 60 J = IA(I), IA(I+1) - 1 DO 60 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
Y(K) = Y(K) + ALPHA*AS(J)*X(I) Y(K) = Y(K) + ALPHA*AS(J)*X(I)
ACC = ACC + AS(J)*X(K) ACC = ACC + AS(J)*X(K)
60 CONTINUE 60 CONTINUE
Y(I) = Y(I) + ALPHA*ACC Y(I) = Y(I) + ALPHA*ACC
80 CONTINUE 80 CONTINUE
ELSE ! Perform computations on conjug(A) ELSE ! Perform computations on conjug(A)
DO 82 I = 1, M DO 82 I = 1, M
ACC = ZERO ACC = ZERO
DO 81 J = IA(I), IA(I+1) - 1 DO 81 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I)
ACC = ACC + CONJG(AS(J)) * X(K) ACC = ACC + CONJG(AS(J)) * X(K)
81 CONTINUE 81 CONTINUE
Y(I) = Y(I) + ALPHA*ACC Y(I) = Y(I) + ALPHA*ACC
82 CONTINUE 82 CONTINUE
ENDIF ENDIF
C C
ELSE IF ( .NOT. UNI) THEN ELSE IF ( .NOT. UNI) THEN
C C
C Check if matrix is lower or upper C Check if matrix is lower or upper
C C
IF ((TRANS.EQ.'L').OR.(TRANS.EQ.'M')) THEN IF ((TRANS.EQ.'L').OR.(TRANS.EQ.'M')) THEN
C C
C LOWER CASE: diagonal element is the last element of row C LOWER CASE: diagonal element is the last element of row
C ....OK! C ....OK!
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 100 I = 1, M DO 100 I = 1, M
Y(I) = BETA*Y(I) Y(I) = BETA*Y(I)
100 CONTINUE 100 CONTINUE
ELSE ELSE
DO I = 1, M DO I = 1, M
Y(I) = ZERO Y(I) = ZERO
ENDDO ENDDO
ENDIF ENDIF
IF (TRANS.EQ.'L') THEN IF (TRANS.EQ.'L') THEN
DO 140 I = 1, M DO 140 I = 1, M
ACC = ZERO ACC = ZERO
DO 120 J = IA(I), IA(I+1) - 1 ! it was -2 DO 120 J = IA(I), IA(I+1) - 1 ! it was -2
K = JA(J) K = JA(J)
C C
C To be optimized C To be optimized
C C
IF (K.NE.I) THEN !for symmetric element IF (K.NE.I) THEN !for symmetric element
Y(K) = Y(K) + ALPHA*AS(J)*X(I) Y(K) = Y(K) + ALPHA*AS(J)*X(I)
ENDIF ENDIF
ACC = ACC + AS(J)*X(K) ACC = ACC + AS(J)*X(K)
120 CONTINUE 120 CONTINUE
Y(I) = Y(I) + ALPHA*ACC Y(I) = Y(I) + ALPHA*ACC
140 CONTINUE 140 CONTINUE
ELSE ! Perform computations on conjug(A) ELSE ! Perform computations on conjug(A)
DO 142 I = 1, M DO 142 I = 1, M
ACC = ZERO ACC = ZERO
DO 141 J = IA(I), IA(I+1) - 1 ! it was -2 DO 141 J = IA(I), IA(I+1) - 1 ! it was -2
K = JA(J) K = JA(J)
C C
C To be optimized C To be optimized
C C
IF (K.NE.I) THEN !for symmetric element IF (K.NE.I) THEN !for symmetric element
Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I)
ENDIF ENDIF
ACC = ACC + CONJG(AS(J)) * X(K) ACC = ACC + CONJG(AS(J)) * X(K)
141 CONTINUE 141 CONTINUE
Y(I) = Y(I) + ALPHA * ACC Y(I) = Y(I) + ALPHA * ACC
142 CONTINUE 142 CONTINUE
ENDIF ENDIF
ELSE ! ....Trans<>L, M ELSE ! ....Trans<>L, M
C C
C UPPER CASE C UPPER CASE
C ....OK!! C ....OK!!
C C
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 160 I = 1, M DO 160 I = 1, M
Y(I) = BETA*Y(I) Y(I) = BETA*Y(I)
160 CONTINUE 160 CONTINUE
ELSE ELSE
DO I = 1, M DO I = 1, M
Y(I) = ZERO Y(I) = ZERO
ENDDO ENDDO
ENDIF ENDIF
IF (TRANS.EQ.'U') THEN IF (TRANS.EQ.'U') THEN
DO 200 I = 1, M DO 200 I = 1, M
ACC = ZERO ACC = ZERO
DO 180 J = IA(I) , IA(I+1) - 1 ! removed +1 DO 180 J = IA(I) , IA(I+1) - 1 ! removed +1
K = JA(J) K = JA(J)
C C
C To be optimized C To be optimized
C C
IF(K.NE.I) THEN IF(K.NE.I) THEN
Y(K) = Y(K) + ALPHA*AS(J)*X(I) Y(K) = Y(K) + ALPHA*AS(J)*X(I)
ENDIF ENDIF
ACC = ACC + AS(J)*X(K) ACC = ACC + AS(J)*X(K)
180 CONTINUE 180 CONTINUE
Y(I) = Y(I) + ALPHA*ACC Y(I) = Y(I) + ALPHA*ACC
200 CONTINUE 200 CONTINUE
ELSE ! Perform computations on conjug(A) ELSE ! Perform computations on conjug(A)
DO 202 I = 1, M DO 202 I = 1, M
ACC = ZERO ACC = ZERO
DO 201 J = IA(I) , IA(I+1) - 1 ! removed +1 DO 201 J = IA(I) , IA(I+1) - 1 ! removed +1
K = JA(J) K = JA(J)
C C
C To be optimized C To be optimized
C C
IF(K.NE.I) THEN IF(K.NE.I) THEN
Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I) Y(K) = Y(K) + ALPHA * CONJG(AS(J)) * X(I)
ENDIF ENDIF
ACC = ACC + CONJG(AS(J)) * X(K) ACC = ACC + CONJG(AS(J)) * X(K)
201 CONTINUE 201 CONTINUE
Y(I) = Y(I) + ALPHA*ACC Y(I) = Y(I) + ALPHA*ACC
202 CONTINUE 202 CONTINUE
ENDIF ENDIF
END IF ! ......TRANS=='L' END IF ! ......TRANS=='L'
END IF ! ......Not UNI END IF ! ......Not UNI
C C
ELSE IF ((.NOT.TRA).AND.(.NOT.COTRA)) THEN !......NOT SYM ELSE IF ((.NOT.TRA).AND.(.NOT.COTRA)) THEN !......NOT SYM
IF ( .NOT. UNI) THEN IF ( .NOT. UNI) THEN
C C
C .......General Not Unit, No Traspose C .......General Not Unit, No Traspose
C C
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 240 I = 1, M DO 240 I = 1, M
ACC = ZERO ACC = ZERO
DO 220 J = IA(I), IA(I+1) - 1 DO 220 J = IA(I), IA(I+1) - 1
ACC = ACC + AS(J)*X(JA(J)) ACC = ACC + AS(J)*X(JA(J))
220 CONTINUE 220 CONTINUE
Y(I) = ALPHA*ACC + BETA*Y(I) Y(I) = ALPHA*ACC + BETA*Y(I)
240 CONTINUE 240 CONTINUE
ELSE ELSE
DO I = 1, M DO I = 1, M
ACC = ZERO ACC = ZERO
DO J = IA(I), IA(I+1) - 1 DO J = IA(I), IA(I+1) - 1
ACC = ACC + AS(J)*X(JA(J)) ACC = ACC + AS(J)*X(JA(J))
ENDDO ENDDO
Y(I) = ALPHA*ACC Y(I) = ALPHA*ACC
ENDDO ENDDO
ENDIF ENDIF
C C
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
C C
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 280 I = 1, M DO 280 I = 1, M
ACC = ZERO ACC = ZERO
DO 260 J = IA(I), IA(I+1) - 1 DO 260 J = IA(I), IA(I+1) - 1
ACC = ACC + AS(J)*X(JA(J)) ACC = ACC + AS(J)*X(JA(J))
260 CONTINUE 260 CONTINUE
Y(I) = ALPHA*(ACC+X(I)) + BETA*Y(I) Y(I) = ALPHA*(ACC+X(I)) + BETA*Y(I)
280 CONTINUE 280 CONTINUE
ELSE !(BETA.EQ.ZERO) ELSE !(BETA.EQ.ZERO)
DO I = 1, M DO I = 1, M
ACC = ZERO ACC = ZERO
DO J = IA(I), IA(I+1) - 1 DO J = IA(I), IA(I+1) - 1
ACC = ACC + AS(J)*X(JA(J)) ACC = ACC + AS(J)*X(JA(J))
ENDDO ENDDO
Y(I) = ALPHA*(ACC+X(I)) Y(I) = ALPHA*(ACC+X(I))
ENDDO ENDDO
ENDIF ENDIF
END IF !....End Testing on UNI END IF !....End Testing on UNI
C C
ELSE IF (TRA) THEN !....Else on SYM (swapped M and N) ELSE IF (TRA) THEN !....Else on SYM (swapped M and N)
C C
IF ( .NOT. UNI) THEN IF ( .NOT. UNI) THEN
C C
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 300 I = 1, M DO 300 I = 1, M
Y(I) = BETA*Y(I) Y(I) = BETA*Y(I)
300 CONTINUE 300 CONTINUE
ELSE !(BETA.EQ.ZERO) ELSE !(BETA.EQ.ZERO)
DO I = 1, M DO I = 1, M
Y(I) = ZERO Y(I) = ZERO
ENDDO ENDDO
ENDIF ENDIF
C C
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
C C
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 320 I = 1, M DO 320 I = 1, M
Y(I) = BETA*Y(I) + ALPHA*X(I) Y(I) = BETA*Y(I) + ALPHA*X(I)
320 CONTINUE 320 CONTINUE
ELSE !(BETA.EQ.ZERO) ELSE !(BETA.EQ.ZERO)
DO I = 1, M DO I = 1, M
Y(I) = ALPHA*X(I) Y(I) = ALPHA*X(I)
ENDDO ENDDO
ENDIF ENDIF
C C
END IF !....UNI END IF !....UNI
C C
IF (ALPHA.EQ.ONE) THEN IF (ALPHA.EQ.ONE) THEN
C C
DO 360 I = 1, N DO 360 I = 1, N
DO 340 J = IA(I), IA(I+1) - 1 DO 340 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
Y(K) = Y(K) + AS(J)*X(I) Y(K) = Y(K) + AS(J)*X(I)
340 CONTINUE 340 CONTINUE
360 CONTINUE 360 CONTINUE
C C
ELSE IF (ALPHA.EQ.-ONE) THEN ELSE IF (ALPHA.EQ.-ONE) THEN
C C
DO 400 I = 1, n DO 400 I = 1, n
DO 380 J = IA(I), IA(I+1) - 1 DO 380 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
Y(K) = Y(K) - AS(J)*X(I) Y(K) = Y(K) - AS(J)*X(I)
380 CONTINUE 380 CONTINUE
400 CONTINUE 400 CONTINUE
C C
ELSE !.....Else on TRA ELSE !.....Else on TRA
C C
C This work array is used for efficiency C This work array is used for efficiency
C C
DO 420 I = 1, N DO 420 I = 1, N
WORK(I) = ALPHA*X(I) WORK(I) = ALPHA*X(I)
420 CONTINUE 420 CONTINUE
C C
DO 460 I = 1, n DO 460 I = 1, n
DO 440 J = IA(I), IA(I+1) - 1 DO 440 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
Y(K) = Y(K) + AS(J)*WORK(I) Y(K) = Y(K) + AS(J)*WORK(I)
440 CONTINUE 440 CONTINUE
460 CONTINUE 460 CONTINUE
C 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 C
IF ( .NOT. UNI) THEN IF ( .NOT. UNI) THEN
C C
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 500 I = 1, M DO 500 I = 1, M
Y(I) = BETA*Y(I) Y(I) = BETA*Y(I)
500 CONTINUE 500 CONTINUE
ELSE !(BETA.EQ.ZERO) ELSE !(BETA.EQ.ZERO)
DO I = 1, M DO I = 1, M
Y(I) = ZERO Y(I) = ZERO
ENDDO ENDDO
ENDIF ENDIF
C C
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
C C
IF (BETA.NE.ZERO) THEN IF (BETA.NE.ZERO) THEN
DO 520 I = 1, M DO 520 I = 1, M
Y(I) = BETA*Y(I) + ALPHA*X(I) Y(I) = BETA*Y(I) + ALPHA*X(I)
520 CONTINUE 520 CONTINUE
ELSE !(BETA.EQ.ZERO) ELSE !(BETA.EQ.ZERO)
DO I = 1, M DO I = 1, M
Y(I) = ALPHA*X(I) Y(I) = ALPHA*X(I)
ENDDO ENDDO
ENDIF ENDIF
C C
END IF !....UNI END IF !....UNI
C C
IF (ALPHA.EQ.ONE) THEN IF (ALPHA.EQ.ONE) THEN
C C
DO 560 I = 1, N DO 560 I = 1, N
DO 540 J = IA(I), IA(I+1) - 1 DO 540 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
Y(K) = Y(K) + CONJG(AS(J))*X(I) Y(K) = Y(K) + CONJG(AS(J))*X(I)
540 CONTINUE 540 CONTINUE
560 CONTINUE 560 CONTINUE
C C
ELSE IF (ALPHA.EQ.-ONE) THEN ELSE IF (ALPHA.EQ.-ONE) THEN
C C
DO 600 I = 1, n DO 600 I = 1, n
DO 580 J = IA(I), IA(I+1) - 1 DO 580 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
Y(K) = Y(K) - CONJG(AS(J))*X(I) Y(K) = Y(K) - CONJG(AS(J))*X(I)
580 CONTINUE 580 CONTINUE
600 CONTINUE 600 CONTINUE
C C
ELSE !.....Else on TRA ELSE !.....Else on TRA
C C
C This work array is used for efficiency C This work array is used for efficiency
C C
DO 620 I = 1, N DO 620 I = 1, N
WORK(I) = ALPHA*X(I) WORK(I) = ALPHA*X(I)
620 CONTINUE 620 CONTINUE
C C
DO 660 I = 1, n DO 660 I = 1, n
DO 640 J = IA(I), IA(I+1) - 1 DO 640 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
Y(K) = Y(K) + CONJG(AS(J))*WORK(I) Y(K) = Y(K) + CONJG(AS(J))*WORK(I)
640 CONTINUE 640 CONTINUE
660 CONTINUE 660 CONTINUE
C 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 C
RETURN RETURN
C C

242
src/serial/csr/zsrsv.f
Unescape Escape View File

@ -13,129 +13,129 @@
COTRA = (TRANS.EQ.'C') COTRA = (TRANS.EQ.'C')
LOW = (UPLO.EQ.'L') LOW = (UPLO.EQ.'L')
IF ((.NOT.TRA).AND.(.NOT.COTRA)) THEN IF ((.NOT.TRA).AND.(.NOT.COTRA)) THEN
IF (LOW) THEN IF (LOW) THEN
IF (.NOT.UNI) THEN IF (.NOT.UNI) THEN
DO 40 I = 1, N DO 40 I = 1, N
ACC = ZERO ACC = ZERO
DO 20 J = IA(I), IA(I+1) - 2 DO 20 J = IA(I), IA(I+1) - 2
ACC = ACC + AS(J)*X(JA(J)) ACC = ACC + AS(J)*X(JA(J))
20 CONTINUE 20 CONTINUE
X(I) = (B(I)-ACC)/AS(IA(I+1)-1) X(I) = (B(I)-ACC)/AS(IA(I+1)-1)
40 CONTINUE 40 CONTINUE
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
DO 80 I = 1, N DO 80 I = 1, N
ACC = ZERO ACC = ZERO
DO 60 J = IA(I), IA(I+1) - 1 DO 60 J = IA(I), IA(I+1) - 1
ACC = ACC + AS(J)*X(JA(J)) ACC = ACC + AS(J)*X(JA(J))
60 CONTINUE 60 CONTINUE
X(I) = B(I) - ACC X(I) = B(I) - ACC
80 CONTINUE 80 CONTINUE
END IF END IF
ELSE IF (.NOT.LOW) THEN ELSE IF (.NOT.LOW) THEN
IF (.NOT.UNI) THEN IF (.NOT.UNI) THEN
DO 120 I = N, 1, -1 DO 120 I = N, 1, -1
ACC = ZERO ACC = ZERO
DO 100 J = IA(I) + 1, IA(I+1) - 1 DO 100 J = IA(I) + 1, IA(I+1) - 1
ACC = ACC + AS(J)*X(JA(J)) ACC = ACC + AS(J)*X(JA(J))
100 CONTINUE 100 CONTINUE
X(I) = (B(I)-ACC)/AS(IA(I)) X(I) = (B(I)-ACC)/AS(IA(I))
120 CONTINUE 120 CONTINUE
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
DO 160 I = N, 1, -1 DO 160 I = N, 1, -1
ACC = ZERO ACC = ZERO
DO 140 J = IA(I), IA(I+1) - 1 DO 140 J = IA(I), IA(I+1) - 1
ACC = ACC + AS(J)*X(JA(J)) ACC = ACC + AS(J)*X(JA(J))
140 CONTINUE 140 CONTINUE
X(I) = B(I) - ACC X(I) = B(I) - ACC
160 CONTINUE 160 CONTINUE
END IF END IF
END IF END IF
ELSE IF (TRA) THEN ELSE IF (TRA) THEN
DO 180 I = 1, N DO 180 I = 1, N
X(I) = B(I) X(I) = B(I)
180 CONTINUE 180 CONTINUE
IF (LOW) THEN IF (LOW) THEN
IF (.NOT.UNI) THEN IF (.NOT.UNI) THEN
DO 220 I = N, 1, -1 DO 220 I = N, 1, -1
X(I) = X(I)/AS(IA(I+1)-1) X(I) = X(I)/AS(IA(I+1)-1)
ACC = X(I) ACC = X(I)
DO 200 J = IA(I), IA(I+1) - 2 DO 200 J = IA(I), IA(I+1) - 2
K = JA(J) K = JA(J)
X(K) = X(K) - AS(J)*ACC X(K) = X(K) - AS(J)*ACC
200 CONTINUE 200 CONTINUE
220 CONTINUE 220 CONTINUE
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
DO 260 I = N, 1, -1 DO 260 I = N, 1, -1
ACC = X(I) ACC = X(I)
DO 240 J = IA(I), IA(I+1) - 1 DO 240 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
X(K) = X(K) - AS(J)*ACC X(K) = X(K) - AS(J)*ACC
240 CONTINUE 240 CONTINUE
260 CONTINUE 260 CONTINUE
END IF END IF
ELSE IF (.NOT.LOW) THEN ELSE IF (.NOT.LOW) THEN
IF (.NOT.UNI) THEN IF (.NOT.UNI) THEN
DO 300 I = 1, N DO 300 I = 1, N
X(I) = X(I)/AS(IA(I)) X(I) = X(I)/AS(IA(I))
ACC = X(I) ACC = X(I)
DO 280 J = IA(I) + 1, IA(I+1) - 1 DO 280 J = IA(I) + 1, IA(I+1) - 1
K = JA(J) K = JA(J)
X(K) = X(K) - AS(J)*ACC X(K) = X(K) - AS(J)*ACC
280 CONTINUE 280 CONTINUE
300 CONTINUE 300 CONTINUE
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
DO 340 I = 1, N DO 340 I = 1, N
ACC = X(I) ACC = X(I)
DO 320 J = IA(I), IA(I+1) - 1 DO 320 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
X(K) = X(K) - AS(J)*ACC X(K) = X(K) - AS(J)*ACC
320 CONTINUE 320 CONTINUE
340 CONTINUE 340 CONTINUE
END IF END IF
END IF END IF
ELSE IF (COTRA) THEN ELSE IF (COTRA) THEN
DO 580 I = 1, N DO 580 I = 1, N
X(I) = B(I) X(I) = B(I)
580 CONTINUE 580 CONTINUE
IF (LOW) THEN IF (LOW) THEN
IF (.NOT.UNI) THEN IF (.NOT.UNI) THEN
DO 620 I = N, 1, -1 DO 620 I = N, 1, -1
X(I) = X(I)/CONJG(AS(IA(I+1)-1)) X(I) = X(I)/CONJG(AS(IA(I+1)-1))
ACC = X(I) ACC = X(I)
DO 600 J = IA(I), IA(I+1) - 2 DO 600 J = IA(I), IA(I+1) - 2
K = JA(J) K = JA(J)
X(K) = X(K) - CONJG(AS(J))*ACC X(K) = X(K) - CONJG(AS(J))*ACC
600 CONTINUE 600 CONTINUE
620 CONTINUE 620 CONTINUE
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
DO 660 I = N, 1, -1 DO 660 I = N, 1, -1
ACC = X(I) ACC = X(I)
DO 640 J = IA(I), IA(I+1) - 1 DO 640 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
X(K) = X(K) - CONJG(AS(J))*ACC X(K) = X(K) - CONJG(AS(J))*ACC
640 CONTINUE 640 CONTINUE
660 CONTINUE 660 CONTINUE
END IF END IF
ELSE IF (.NOT.LOW) THEN ELSE IF (.NOT.LOW) THEN
IF (.NOT.UNI) THEN IF (.NOT.UNI) THEN
DO 700 I = 1, N DO 700 I = 1, N
X(I) = X(I)/CONJG(AS(IA(I))) X(I) = X(I)/CONJG(AS(IA(I)))
ACC = X(I) ACC = X(I)
DO 680 J = IA(I) + 1, IA(I+1) - 1 DO 680 J = IA(I) + 1, IA(I+1) - 1
K = JA(J) K = JA(J)
X(K) = X(K) - CONJG(AS(J))*ACC X(K) = X(K) - CONJG(AS(J))*ACC
680 CONTINUE 680 CONTINUE
700 CONTINUE 700 CONTINUE
ELSE IF (UNI) THEN ELSE IF (UNI) THEN
DO 740 I = 1, N DO 740 I = 1, N
ACC = X(I) ACC = X(I)
DO 720 J = IA(I), IA(I+1) - 1 DO 720 J = IA(I), IA(I+1) - 1
K = JA(J) K = JA(J)
X(K) = X(K) - CONJG(AS(J))*ACC X(K) = X(K) - CONJG(AS(J))*ACC
720 CONTINUE 720 CONTINUE
740 CONTINUE 740 CONTINUE
END IF END IF
END IF END IF
END IF END IF
RETURN RETURN
END END

2
src/serial/f77/zcssm.f
Unescape Escape View File

@ -231,6 +231,7 @@ C
C Error handling C Error handling
C C
IF(IERROR.NE.0) THEN IF(IERROR.NE.0) THEN
write(0,*) 'zcssm ierror',ierror
CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL) CALL FCPSB_ERRPUSH(IERROR,NAME,INT_VAL)
GOTO 9999 GOTO 9999
@ -335,6 +336,7 @@ C
RETURN RETURN
9999 CONTINUE 9999 CONTINUE
CALL FCPSB_ERRACTIONRESTORE(ERR_ACT) CALL FCPSB_ERRACTIONRESTORE(ERR_ACT)
IF ( ERR_ACT .NE. 0 ) THEN IF ( ERR_ACT .NE. 0 ) THEN
CALL FCPSB_SERROR() CALL FCPSB_SERROR()

2
src/serial/psb_zcssv.f90
Unescape Escape View File

@ -80,7 +80,7 @@ subroutine psb_zcssv(alpha,t,b,beta,c,info,trans,unitd,d)
call zcssm(lt,m,n,alpha,lu,ddl,& call zcssm(lt,m,n,alpha,lu,ddl,&
& t%pl,t%fida,t%descra,t%aspk,t%ia1,t%ia2,t%infoa,t%pr,& & t%pl,t%fida,t%descra,t%aspk,t%ia1,t%ia2,t%infoa,t%pr,&
& b,lb,beta,c,lc,work,iwsz,info) & b,lb,beta,c,lc,work,iwsz,info)
if (.not.present(d)) then if (.not.present(d)) then
deallocate(ddl) deallocate(ddl)
endif endif

Write Preview
Loading…
Cancel
Save