First batch of Fortran to C conversions
This commit is contained in:
Adam Parler
parent
f33185ce7f
commit
afee077224
29 changed files with 2566 additions and 2674 deletions
2
.gitignore
vendored
Normal file
2
.gitignore
vendored
Normal file
|
|
@ -0,0 +1,2 @@
|
|||
.vscode/*
|
||||
build/*
|
||||
33
blas.h
Normal file
33
blas.h
Normal file
|
|
@ -0,0 +1,33 @@
|
|||
#ifndef _BLAS_H
|
||||
#define _BLAS_H
|
||||
|
||||
#ifdef __cplusplus
|
||||
extern "C" {
|
||||
#endif
|
||||
|
||||
#include <complex.h>
|
||||
|
||||
void caxpy ( int n, float complex ca, float complex cx[], int incx, float complex cy[], int incy );
|
||||
void ccopy ( int n, float complex cx[], int incx, float complex cy[], int incy );
|
||||
float complex cdotc ( int n, float complex cx[], int incx, float complex cy[], int incy );
|
||||
float complex cdotu ( int n, float complex cx[], int incx, float complex cy[], int incy );
|
||||
void cgemm(char transa, char transb, int m, int n, int k, complex float alpha,
|
||||
complex float a[], int lda, complex float b[], int ldb, complex float beta,
|
||||
complex float c[], int ldc);
|
||||
void crotg ( float complex *ca, float complex cb, float *c, float complex *s );
|
||||
void cscal ( int n, float complex ca, float complex cx[], int incx );
|
||||
void csrot ( int n, float complex cx[], int incx, float complex cy[], int incy, float c, float s );
|
||||
void csscal ( int n, float sa, float complex cx[], int incx );
|
||||
void cswap ( int n, float complex cx[], int incx, float complex cy[], int incy );
|
||||
int icamax ( int n, float complex x[], int incx );
|
||||
float scasum ( int n, float complex x[], int incx );
|
||||
float scnrm2 ( int n, float complex x[], int incx );
|
||||
|
||||
|
||||
|
||||
#ifdef __cplusplus
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
#endif
|
||||
106
src/caxpy.c
106
src/caxpy.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CAXPY
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -84,56 +89,51 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CAXPY(N,CA,CX,INCX,CY,INCY)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
COMPLEX CA
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*),CY(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
INTEGER I,IX,IY
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
REAL SCABS1
|
||||
EXTERNAL SCABS1
|
||||
* ..
|
||||
IF (N.LE.0) RETURN
|
||||
IF (SCABS1(CA).EQ.0.0E+0) RETURN
|
||||
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
|
||||
*
|
||||
* code for both increments equal to 1
|
||||
*
|
||||
DO I = 1,N
|
||||
CY(I) = CY(I) + CA*CX(I)
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for unequal increments or equal increments
|
||||
* not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
IY = 1
|
||||
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
|
||||
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
|
||||
DO I = 1,N
|
||||
CY(IY) = CY(IY) + CA*CX(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of CAXPY
|
||||
*
|
||||
END
|
||||
*/
|
||||
void caxpy ( int n, float complex ca, float complex cx[], int incx,
|
||||
float complex cy[], int incy ) {
|
||||
//
|
||||
// -- Reference BLAS level1 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
//
|
||||
int i, ix, iy;
|
||||
|
||||
if ( n <= 0 ) return;
|
||||
|
||||
|
||||
if ( scabs1 ( ca ) == 0.0 ) return;
|
||||
|
||||
if ( incx != 1 || incy != 1 ) {
|
||||
if ( 0 <= incx ) {
|
||||
ix = 0;
|
||||
}
|
||||
else {
|
||||
ix = ( -n + 1 ) * incx;
|
||||
}
|
||||
|
||||
if ( 0 <= incy ) {
|
||||
iy = 0;
|
||||
}
|
||||
else {
|
||||
iy = ( -n + 1 ) * incy;
|
||||
}
|
||||
|
||||
for (i = 0; i < n; i++) {
|
||||
cy[iy] = cy[iy] + ca * cx[ix];
|
||||
ix = ix + incx;
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
else {
|
||||
for (i = 0; i < n; i++) {
|
||||
cy[i] = cy[i] + ca * cx[i];
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of CAXPY
|
||||
|
||||
}
|
||||
|
|
|
|||
94
src/ccopy.c
94
src/ccopy.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CCOPY
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -77,49 +82,46 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CCOPY(N,CX,INCX,CY,INCY)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*),CY(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
INTEGER I,IX,IY
|
||||
* ..
|
||||
IF (N.LE.0) RETURN
|
||||
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
|
||||
*
|
||||
* code for both increments equal to 1
|
||||
*
|
||||
DO I = 1,N
|
||||
CY(I) = CX(I)
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for unequal increments or equal increments
|
||||
* not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
IY = 1
|
||||
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
|
||||
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
|
||||
DO I = 1,N
|
||||
CY(IY) = CX(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of CCOPY
|
||||
*
|
||||
END
|
||||
*/
|
||||
void ccopy ( int n, float complex cx[], int incx, float complex cy[],
|
||||
int incy ) {
|
||||
|
||||
// -- Reference BLAS level1 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
|
||||
int i, ix, iy;
|
||||
|
||||
if ( n <= 0 ) {
|
||||
return;
|
||||
}
|
||||
|
||||
if ( incx != 1 || incy != 1 ) {
|
||||
if ( 0 <= incx ) {
|
||||
ix = 0;
|
||||
} else {
|
||||
ix = ( -n + 1 ) * incx;
|
||||
}
|
||||
|
||||
if ( 0 <= incy ) {
|
||||
iy = 0;
|
||||
} else {
|
||||
iy = ( -n + 1 ) * incy;
|
||||
}
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
cy[iy] = cx[ix];
|
||||
ix = ix + incx;
|
||||
iy = iy + incy;
|
||||
}
|
||||
} else {
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
cy[i] = cx[i];
|
||||
}
|
||||
}
|
||||
return;
|
||||
|
||||
// End of CCOPY
|
||||
|
||||
}
|
||||
|
|
|
|||
106
src/cdotc.c
106
src/cdotc.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CDOTC
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -79,56 +84,51 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*),CY(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
COMPLEX CTEMP
|
||||
INTEGER I,IX,IY
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC CONJG
|
||||
* ..
|
||||
CTEMP = (0.0,0.0)
|
||||
CDOTC = (0.0,0.0)
|
||||
IF (N.LE.0) RETURN
|
||||
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
|
||||
*
|
||||
* code for both increments equal to 1
|
||||
*
|
||||
DO I = 1,N
|
||||
CTEMP = CTEMP + CONJG(CX(I))*CY(I)
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for unequal increments or equal increments
|
||||
* not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
IY = 1
|
||||
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
|
||||
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
|
||||
DO I = 1,N
|
||||
CTEMP = CTEMP + CONJG(CX(IX))*CY(IY)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
END DO
|
||||
END IF
|
||||
CDOTC = CTEMP
|
||||
RETURN
|
||||
*
|
||||
* End of CDOTC
|
||||
*
|
||||
END
|
||||
*/
|
||||
float complex cdotc ( int n, float complex cx[], int incx,
|
||||
float complex cy[], int incy ) {
|
||||
|
||||
// -- Reference BLAS level1 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
|
||||
|
||||
int i, ix, iy;
|
||||
float complex value;
|
||||
|
||||
value = 0.0;
|
||||
|
||||
if ( n <= 0 ) {
|
||||
return value;
|
||||
}
|
||||
|
||||
if ( incx == 1 && incy == 1 ) {
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
value = value + ( ~cx[i] ) * cy[i];
|
||||
}
|
||||
} else {
|
||||
if ( 0 <= incx ) {
|
||||
ix = 0;
|
||||
} else {
|
||||
ix = ( -n + 1 ) * incx;
|
||||
}
|
||||
|
||||
if ( 0 <= incy ) {
|
||||
iy = 0;
|
||||
} else {
|
||||
iy = ( -n + 1 ) * incy;
|
||||
}
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
value = value + ( ~cx[ix] ) * cy[iy];
|
||||
ix = ix + incx;
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
|
||||
return value;
|
||||
|
||||
// End of CDOTC
|
||||
|
||||
}
|
||||
|
|
|
|||
103
src/cdotu.c
103
src/cdotu.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CDOTU
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -79,53 +84,51 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
COMPLEX FUNCTION CDOTU(N,CX,INCX,CY,INCY)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*),CY(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
COMPLEX CTEMP
|
||||
INTEGER I,IX,IY
|
||||
* ..
|
||||
CTEMP = (0.0,0.0)
|
||||
CDOTU = (0.0,0.0)
|
||||
IF (N.LE.0) RETURN
|
||||
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
|
||||
*
|
||||
* code for both increments equal to 1
|
||||
*
|
||||
DO I = 1,N
|
||||
CTEMP = CTEMP + CX(I)*CY(I)
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for unequal increments or equal increments
|
||||
* not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
IY = 1
|
||||
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
|
||||
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
|
||||
DO I = 1,N
|
||||
CTEMP = CTEMP + CX(IX)*CY(IY)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
END DO
|
||||
END IF
|
||||
CDOTU = CTEMP
|
||||
RETURN
|
||||
*
|
||||
* End of CDOTU
|
||||
*
|
||||
END
|
||||
*/
|
||||
float complex cdotu ( int n, float complex cx[], int incx,
|
||||
float complex cy[], int incy ) {
|
||||
|
||||
// -- Reference BLAS level1 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
|
||||
|
||||
int i, ix, iy;
|
||||
float complex value;
|
||||
|
||||
value = 0.0;
|
||||
|
||||
if ( n <= 0 ) {
|
||||
return value;
|
||||
}
|
||||
|
||||
if ( incx == 1 && incy == 1 ) {
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
value = value + cx[i] * cy[i];
|
||||
}
|
||||
} else {
|
||||
if ( 0 <= incx ) {
|
||||
ix = 0;
|
||||
} else {
|
||||
ix = ( -n + 1 ) * incx;
|
||||
}
|
||||
|
||||
if ( 0 <= incy ) {
|
||||
iy = 0;
|
||||
} else {
|
||||
iy = ( -n + 1 ) * incy;
|
||||
}
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
value = value + cx[ix] * cy[iy];
|
||||
ix = ix + incx;
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
|
||||
return value;
|
||||
|
||||
// End of CDOTU
|
||||
|
||||
}
|
||||
|
|
|
|||
616
src/cgemm.c
616
src/cgemm.c
|
|
@ -1,3 +1,11 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <complex.h>
|
||||
#include <stdbool.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CGEMM
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -183,295 +191,319 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
|
||||
*
|
||||
* -- Reference BLAS level3 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
COMPLEX ALPHA,BETA
|
||||
INTEGER K,LDA,LDB,LDC,M,N
|
||||
CHARACTER TRANSA,TRANSB
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC CONJG,MAX
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
COMPLEX TEMP
|
||||
INTEGER I,INFO,J,L,NROWA,NROWB
|
||||
LOGICAL CONJA,CONJB,NOTA,NOTB
|
||||
* ..
|
||||
* .. Parameters ..
|
||||
COMPLEX ONE
|
||||
PARAMETER (ONE= (1.0E+0,0.0E+0))
|
||||
COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0E+0,0.0E+0))
|
||||
* ..
|
||||
*
|
||||
* Set NOTA and NOTB as true if A and B respectively are not
|
||||
* conjugated or transposed, set CONJA and CONJB as true if A and
|
||||
* B respectively are to be transposed but not conjugated and set
|
||||
* NROWA and NROWB as the number of rows of A and B respectively.
|
||||
*
|
||||
NOTA = LSAME(TRANSA,'N')
|
||||
NOTB = LSAME(TRANSB,'N')
|
||||
CONJA = LSAME(TRANSA,'C')
|
||||
CONJB = LSAME(TRANSB,'C')
|
||||
IF (NOTA) THEN
|
||||
NROWA = M
|
||||
ELSE
|
||||
NROWA = K
|
||||
END IF
|
||||
IF (NOTB) THEN
|
||||
NROWB = K
|
||||
ELSE
|
||||
NROWB = N
|
||||
END IF
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
IF ((.NOT.NOTA) .AND. (.NOT.CONJA) .AND.
|
||||
+ (.NOT.LSAME(TRANSA,'T'))) THEN
|
||||
INFO = 1
|
||||
ELSE IF ((.NOT.NOTB) .AND. (.NOT.CONJB) .AND.
|
||||
+ (.NOT.LSAME(TRANSB,'T'))) THEN
|
||||
INFO = 2
|
||||
ELSE IF (M.LT.0) THEN
|
||||
INFO = 3
|
||||
ELSE IF (N.LT.0) THEN
|
||||
INFO = 4
|
||||
ELSE IF (K.LT.0) THEN
|
||||
INFO = 5
|
||||
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
|
||||
INFO = 8
|
||||
ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
|
||||
INFO = 10
|
||||
ELSE IF (LDC.LT.MAX(1,M)) THEN
|
||||
INFO = 13
|
||||
END IF
|
||||
IF (INFO.NE.0) THEN
|
||||
CALL XERBLA('CGEMM ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
|
||||
+ (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
|
||||
*
|
||||
* And when alpha.eq.zero.
|
||||
*
|
||||
IF (ALPHA.EQ.ZERO) THEN
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
DO 20 J = 1,N
|
||||
DO 10 I = 1,M
|
||||
C(I,J) = ZERO
|
||||
10 CONTINUE
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
DO 40 J = 1,N
|
||||
DO 30 I = 1,M
|
||||
C(I,J) = BETA*C(I,J)
|
||||
30 CONTINUE
|
||||
40 CONTINUE
|
||||
END IF
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Start the operations.
|
||||
*
|
||||
IF (NOTB) THEN
|
||||
IF (NOTA) THEN
|
||||
*
|
||||
* Form C := alpha*A*B + beta*C.
|
||||
*
|
||||
DO 90 J = 1,N
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
DO 50 I = 1,M
|
||||
C(I,J) = ZERO
|
||||
50 CONTINUE
|
||||
ELSE IF (BETA.NE.ONE) THEN
|
||||
DO 60 I = 1,M
|
||||
C(I,J) = BETA*C(I,J)
|
||||
60 CONTINUE
|
||||
END IF
|
||||
DO 80 L = 1,K
|
||||
TEMP = ALPHA*B(L,J)
|
||||
DO 70 I = 1,M
|
||||
C(I,J) = C(I,J) + TEMP*A(I,L)
|
||||
70 CONTINUE
|
||||
80 CONTINUE
|
||||
90 CONTINUE
|
||||
ELSE IF (CONJA) THEN
|
||||
*
|
||||
* Form C := alpha*A**H*B + beta*C.
|
||||
*
|
||||
DO 120 J = 1,N
|
||||
DO 110 I = 1,M
|
||||
TEMP = ZERO
|
||||
DO 100 L = 1,K
|
||||
TEMP = TEMP + CONJG(A(L,I))*B(L,J)
|
||||
100 CONTINUE
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
C(I,J) = ALPHA*TEMP
|
||||
ELSE
|
||||
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
|
||||
END IF
|
||||
110 CONTINUE
|
||||
120 CONTINUE
|
||||
ELSE
|
||||
*
|
||||
* Form C := alpha*A**T*B + beta*C
|
||||
*
|
||||
DO 150 J = 1,N
|
||||
DO 140 I = 1,M
|
||||
TEMP = ZERO
|
||||
DO 130 L = 1,K
|
||||
TEMP = TEMP + A(L,I)*B(L,J)
|
||||
130 CONTINUE
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
C(I,J) = ALPHA*TEMP
|
||||
ELSE
|
||||
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
|
||||
END IF
|
||||
140 CONTINUE
|
||||
150 CONTINUE
|
||||
END IF
|
||||
ELSE IF (NOTA) THEN
|
||||
IF (CONJB) THEN
|
||||
*
|
||||
* Form C := alpha*A*B**H + beta*C.
|
||||
*
|
||||
DO 200 J = 1,N
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
DO 160 I = 1,M
|
||||
C(I,J) = ZERO
|
||||
160 CONTINUE
|
||||
ELSE IF (BETA.NE.ONE) THEN
|
||||
DO 170 I = 1,M
|
||||
C(I,J) = BETA*C(I,J)
|
||||
170 CONTINUE
|
||||
END IF
|
||||
DO 190 L = 1,K
|
||||
TEMP = ALPHA*CONJG(B(J,L))
|
||||
DO 180 I = 1,M
|
||||
C(I,J) = C(I,J) + TEMP*A(I,L)
|
||||
180 CONTINUE
|
||||
190 CONTINUE
|
||||
200 CONTINUE
|
||||
ELSE
|
||||
*
|
||||
* Form C := alpha*A*B**T + beta*C
|
||||
*
|
||||
DO 250 J = 1,N
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
DO 210 I = 1,M
|
||||
C(I,J) = ZERO
|
||||
210 CONTINUE
|
||||
ELSE IF (BETA.NE.ONE) THEN
|
||||
DO 220 I = 1,M
|
||||
C(I,J) = BETA*C(I,J)
|
||||
220 CONTINUE
|
||||
END IF
|
||||
DO 240 L = 1,K
|
||||
TEMP = ALPHA*B(J,L)
|
||||
DO 230 I = 1,M
|
||||
C(I,J) = C(I,J) + TEMP*A(I,L)
|
||||
230 CONTINUE
|
||||
240 CONTINUE
|
||||
250 CONTINUE
|
||||
END IF
|
||||
ELSE IF (CONJA) THEN
|
||||
IF (CONJB) THEN
|
||||
*
|
||||
* Form C := alpha*A**H*B**H + beta*C.
|
||||
*
|
||||
DO 280 J = 1,N
|
||||
DO 270 I = 1,M
|
||||
TEMP = ZERO
|
||||
DO 260 L = 1,K
|
||||
TEMP = TEMP + CONJG(A(L,I))*CONJG(B(J,L))
|
||||
260 CONTINUE
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
C(I,J) = ALPHA*TEMP
|
||||
ELSE
|
||||
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
|
||||
END IF
|
||||
270 CONTINUE
|
||||
280 CONTINUE
|
||||
ELSE
|
||||
*
|
||||
* Form C := alpha*A**H*B**T + beta*C
|
||||
*
|
||||
DO 310 J = 1,N
|
||||
DO 300 I = 1,M
|
||||
TEMP = ZERO
|
||||
DO 290 L = 1,K
|
||||
TEMP = TEMP + CONJG(A(L,I))*B(J,L)
|
||||
290 CONTINUE
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
C(I,J) = ALPHA*TEMP
|
||||
ELSE
|
||||
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
|
||||
END IF
|
||||
300 CONTINUE
|
||||
310 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IF (CONJB) THEN
|
||||
*
|
||||
* Form C := alpha*A**T*B**H + beta*C
|
||||
*
|
||||
DO 340 J = 1,N
|
||||
DO 330 I = 1,M
|
||||
TEMP = ZERO
|
||||
DO 320 L = 1,K
|
||||
TEMP = TEMP + A(L,I)*CONJG(B(J,L))
|
||||
320 CONTINUE
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
C(I,J) = ALPHA*TEMP
|
||||
ELSE
|
||||
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
|
||||
END IF
|
||||
330 CONTINUE
|
||||
340 CONTINUE
|
||||
ELSE
|
||||
*
|
||||
* Form C := alpha*A**T*B**T + beta*C
|
||||
*
|
||||
DO 370 J = 1,N
|
||||
DO 360 I = 1,M
|
||||
TEMP = ZERO
|
||||
DO 350 L = 1,K
|
||||
TEMP = TEMP + A(L,I)*B(J,L)
|
||||
350 CONTINUE
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
C(I,J) = ALPHA*TEMP
|
||||
ELSE
|
||||
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
|
||||
END IF
|
||||
360 CONTINUE
|
||||
370 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of CGEMM
|
||||
*
|
||||
END
|
||||
*/
|
||||
void cgemm(char transa, char transb, int m, int n, int k, complex float alpha,
|
||||
complex float a[], int lda, complex float b[], int ldb, complex float beta,
|
||||
complex float c[], int ldc) {
|
||||
|
||||
// -- Reference BLAS level3 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
|
||||
// Local Scalars ..
|
||||
int i, j, l, nrowa, nrowb, info, ncola;
|
||||
complex float temp;
|
||||
bool conja, conjb, nota, notb;
|
||||
|
||||
complex float one = (1.0e0, 0.0e0);
|
||||
complex float zero = (0.0e0, 0.0e0);
|
||||
|
||||
// Set NOTA and NOTB as true if A and B respectively are not
|
||||
// conjugated or transposed, set CONJA and CONJB as true if A and
|
||||
// B respectively are to be transposed but not conjugated and set
|
||||
// NROWA, NCOLA and NROWB as the number of rows and columns of A
|
||||
// and the number of rows of B respectively.
|
||||
|
||||
nota = ( ( transa == 'N' ) || ( transa == 'n' ) );
|
||||
notb = ( ( transb == 'N' ) || ( transb == 'n' ) );
|
||||
conja = ( ( transa == 'C' ) || ( transa == 'c' ) );
|
||||
conjb = ( ( transb == 'C' ) || ( transb == 'c' ) );
|
||||
|
||||
if ( nota ) {
|
||||
nrowa = m;
|
||||
ncola = k;
|
||||
}
|
||||
else {
|
||||
nrowa = k;
|
||||
ncola = m;
|
||||
}
|
||||
|
||||
if (notb) {
|
||||
nrowb = k;
|
||||
}
|
||||
else {
|
||||
nrowb = n;
|
||||
}
|
||||
|
||||
// Test the input parameters.
|
||||
|
||||
info = 0;
|
||||
if ( ! ( transa == 'N' || transa == 'n' ||
|
||||
transa == 'C' || transa == 'c' ||
|
||||
transa == 'T' || transa == 't' ) ) {
|
||||
info = 1;
|
||||
}
|
||||
|
||||
if ( ! ( transb == 'N' || transb == 'n' ||
|
||||
transb == 'C' || transb == 'c' ||
|
||||
transb == 'T' || transb == 't' ) ) {
|
||||
info = 2;
|
||||
}
|
||||
|
||||
if ( m < 0 ) {
|
||||
info = 3;
|
||||
}
|
||||
else if ( n < 0 ) {
|
||||
info = 4;
|
||||
}
|
||||
else if ( k < 0 ) {
|
||||
info = 5;
|
||||
}
|
||||
else if ( lda < i4_max( 1, nrowa ) ) {
|
||||
info = 8;
|
||||
}
|
||||
else if ( ldb < i4_max( 1, nrowb ) ) {
|
||||
info = 10;
|
||||
}
|
||||
else if ( ldc < i4_max( 1, m ) ) {
|
||||
info = 13;
|
||||
}
|
||||
|
||||
if (info != 0) {
|
||||
xerbla('cgemm ', info);
|
||||
return;
|
||||
}
|
||||
|
||||
// Quick return if possible.
|
||||
|
||||
if (( m == 0 ) || (n == 0) || (((alpha == zero) || ( k == 0 )) && ( beta == one ) ) ) {
|
||||
return;
|
||||
}
|
||||
|
||||
// And when alpha.eq.zero.
|
||||
|
||||
if ( alpha == zero ) {
|
||||
if ( beta == zero ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = zero;
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = beta * c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
// Start the operations.
|
||||
|
||||
if ( notb ) {
|
||||
if ( nota ) {
|
||||
|
||||
// Form C := alpha*A*B + beta*C.
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( beta == zero ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = zero;
|
||||
}
|
||||
}
|
||||
else if ( beta != one ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
if ( b[l+j*ldb] != zero ) {
|
||||
temp = alpha * b[l+j*ldb];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = c[i+j*ldc] + temp * a[i+l*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (conja) {
|
||||
|
||||
// Form C := alpha*A**H*B + beta*C.
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = zero;
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
temp = temp + (~a[l+i*lda])* b[l+j*ldb];
|
||||
}
|
||||
if (beta == zero) {
|
||||
c[i+j*ldc] = alpha*temp;
|
||||
}
|
||||
else {
|
||||
c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
|
||||
// Form C := alpha*A**T*B + beta*C
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = zero;
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
temp = temp + a[l+i*lda] * b[l+j*ldb];
|
||||
}
|
||||
if (beta == zero) {
|
||||
c[i+j*ldc] = alpha*temp;
|
||||
}
|
||||
else {
|
||||
c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (nota) {
|
||||
if (conjb) {
|
||||
|
||||
// Form C := alpha*A*B**H + beta*C.
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if (beta == zero) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = zero;
|
||||
}
|
||||
}
|
||||
else if (beta != one) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
if (b[j+l*ldb] != zero) {
|
||||
temp = alpha*(~b[j+l*ldb]);
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = c[i+j*ldc] + temp*a[i+l*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
|
||||
// Form C := alpha*A*B**T + beta*C
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if (beta == zero) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = zero;
|
||||
}
|
||||
}
|
||||
else if (beta != one) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
if (b[j+l*ldb] != zero) {
|
||||
temp = alpha*b[j+l*ldb];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
c[i+j*ldc] = c[i+j*ldc] + temp*a[i+l*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else if (conja) {
|
||||
if (conjb) {
|
||||
|
||||
// Form C := alpha*A**H*B**H + beta*C.
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = zero;
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
temp = temp + (~a[l+i*lda])*(~b[j+l*ldb]);
|
||||
}
|
||||
if (beta == zero) {
|
||||
c[i+j*ldc] = alpha*temp;
|
||||
}
|
||||
else {
|
||||
c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
|
||||
// Form C := alpha*A**H*B**T + beta*C
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = zero;
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
temp = temp + (~a[l+i*lda])*b[j+l*ldb];
|
||||
}
|
||||
if (beta == zero) {
|
||||
c[i+j*ldc] = alpha*temp;
|
||||
}
|
||||
else {
|
||||
c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
if (conjb) {
|
||||
|
||||
// Form C := alpha*A**T*B**H + beta*C
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = zero;
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
temp = temp + a[l+i*lda]*(~b[j+l*ldb]);
|
||||
}
|
||||
if (beta == zero) {
|
||||
c[i+j*ldc] = alpha*temp;
|
||||
}
|
||||
else {
|
||||
c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
|
||||
// Form C := alpha*A**T*B**T + beta*C
|
||||
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = zero;
|
||||
for ( l = 0; l < k; l++ ) {
|
||||
temp = temp + a[l+i*lda]*b[j+l*ldb];
|
||||
}
|
||||
if (beta == zero) {
|
||||
c[i+j*ldc] = alpha*temp;
|
||||
}
|
||||
else {
|
||||
c[i+j*ldc] = alpha*temp + beta*c[i+j*ldc];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of CGEMM .
|
||||
|
||||
}
|
||||
379
src/cgemv.c
379
src/cgemv.c
|
|
@ -1,3 +1,9 @@
|
|||
#include <complex.h>
|
||||
#include <stdbool.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CGEMV
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -154,194 +160,185 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
|
||||
*
|
||||
* -- Reference BLAS level2 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
COMPLEX ALPHA,BETA
|
||||
INTEGER INCX,INCY,LDA,M,N
|
||||
CHARACTER TRANS
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX A(LDA,*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
COMPLEX ONE
|
||||
PARAMETER (ONE= (1.0E+0,0.0E+0))
|
||||
COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0E+0,0.0E+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
COMPLEX TEMP
|
||||
INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY
|
||||
LOGICAL NOCONJ
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC CONJG,MAX
|
||||
* ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
|
||||
+ .NOT.LSAME(TRANS,'C')) THEN
|
||||
INFO = 1
|
||||
ELSE IF (M.LT.0) THEN
|
||||
INFO = 2
|
||||
ELSE IF (N.LT.0) THEN
|
||||
INFO = 3
|
||||
ELSE IF (LDA.LT.MAX(1,M)) THEN
|
||||
INFO = 6
|
||||
ELSE IF (INCX.EQ.0) THEN
|
||||
INFO = 8
|
||||
ELSE IF (INCY.EQ.0) THEN
|
||||
INFO = 11
|
||||
END IF
|
||||
IF (INFO.NE.0) THEN
|
||||
CALL XERBLA('CGEMV ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
|
||||
+ ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
|
||||
*
|
||||
NOCONJ = LSAME(TRANS,'T')
|
||||
*
|
||||
* Set LENX and LENY, the lengths of the vectors x and y, and set
|
||||
* up the start points in X and Y.
|
||||
*
|
||||
IF (LSAME(TRANS,'N')) THEN
|
||||
LENX = N
|
||||
LENY = M
|
||||
ELSE
|
||||
LENX = M
|
||||
LENY = N
|
||||
END IF
|
||||
IF (INCX.GT.0) THEN
|
||||
KX = 1
|
||||
ELSE
|
||||
KX = 1 - (LENX-1)*INCX
|
||||
END IF
|
||||
IF (INCY.GT.0) THEN
|
||||
KY = 1
|
||||
ELSE
|
||||
KY = 1 - (LENY-1)*INCY
|
||||
END IF
|
||||
*
|
||||
* Start the operations. In this version the elements of A are
|
||||
* accessed sequentially with one pass through A.
|
||||
*
|
||||
* First form y := beta*y.
|
||||
*
|
||||
IF (BETA.NE.ONE) THEN
|
||||
IF (INCY.EQ.1) THEN
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
DO 10 I = 1,LENY
|
||||
Y(I) = ZERO
|
||||
10 CONTINUE
|
||||
ELSE
|
||||
DO 20 I = 1,LENY
|
||||
Y(I) = BETA*Y(I)
|
||||
20 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IY = KY
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
DO 30 I = 1,LENY
|
||||
Y(IY) = ZERO
|
||||
IY = IY + INCY
|
||||
30 CONTINUE
|
||||
ELSE
|
||||
DO 40 I = 1,LENY
|
||||
Y(IY) = BETA*Y(IY)
|
||||
IY = IY + INCY
|
||||
40 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
IF (ALPHA.EQ.ZERO) RETURN
|
||||
IF (LSAME(TRANS,'N')) THEN
|
||||
*
|
||||
* Form y := alpha*A*x + y.
|
||||
*
|
||||
JX = KX
|
||||
IF (INCY.EQ.1) THEN
|
||||
DO 60 J = 1,N
|
||||
TEMP = ALPHA*X(JX)
|
||||
DO 50 I = 1,M
|
||||
Y(I) = Y(I) + TEMP*A(I,J)
|
||||
50 CONTINUE
|
||||
JX = JX + INCX
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
DO 80 J = 1,N
|
||||
TEMP = ALPHA*X(JX)
|
||||
IY = KY
|
||||
DO 70 I = 1,M
|
||||
Y(IY) = Y(IY) + TEMP*A(I,J)
|
||||
IY = IY + INCY
|
||||
70 CONTINUE
|
||||
JX = JX + INCX
|
||||
80 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
*
|
||||
* Form y := alpha*A**T*x + y or y := alpha*A**H*x + y.
|
||||
*
|
||||
JY = KY
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 110 J = 1,N
|
||||
TEMP = ZERO
|
||||
IF (NOCONJ) THEN
|
||||
DO 90 I = 1,M
|
||||
TEMP = TEMP + A(I,J)*X(I)
|
||||
90 CONTINUE
|
||||
ELSE
|
||||
DO 100 I = 1,M
|
||||
TEMP = TEMP + CONJG(A(I,J))*X(I)
|
||||
100 CONTINUE
|
||||
END IF
|
||||
Y(JY) = Y(JY) + ALPHA*TEMP
|
||||
JY = JY + INCY
|
||||
110 CONTINUE
|
||||
ELSE
|
||||
DO 140 J = 1,N
|
||||
TEMP = ZERO
|
||||
IX = KX
|
||||
IF (NOCONJ) THEN
|
||||
DO 120 I = 1,M
|
||||
TEMP = TEMP + A(I,J)*X(IX)
|
||||
IX = IX + INCX
|
||||
120 CONTINUE
|
||||
ELSE
|
||||
DO 130 I = 1,M
|
||||
TEMP = TEMP + CONJG(A(I,J))*X(IX)
|
||||
IX = IX + INCX
|
||||
130 CONTINUE
|
||||
END IF
|
||||
Y(JY) = Y(JY) + ALPHA*TEMP
|
||||
JY = JY + INCY
|
||||
140 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of CGEMV
|
||||
*
|
||||
END
|
||||
*/
|
||||
void cgemv(char trans, int m, int n, complex float alpha, complex float a[], int lda, complex float x[], int incx, complex float beta, complex float y[], int incy) {
|
||||
|
||||
// -- Reference BLAS level2 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
|
||||
|
||||
// ..
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Parameters ..
|
||||
complex float one = (1.0E+0,0.0E+0 * I);
|
||||
complex float zero = (0.0E+0,0.0E+0 * I);
|
||||
// ..
|
||||
// .. Local Scalars ..
|
||||
complex float temp;
|
||||
int i,info,ix,iy,j,jx,jy,kx,ky,lenx,leny;
|
||||
bool noconj;
|
||||
|
||||
// ..
|
||||
|
||||
// Test the input parameters.
|
||||
|
||||
info = 0;
|
||||
if (!lsame(trans,'N') && !lsame(trans,'T') &&
|
||||
+ !lsame(trans,'C')) {
|
||||
info = 1;
|
||||
} else if (m < 0) {
|
||||
info = 2;
|
||||
} else if (n < 0) {
|
||||
info = 3;
|
||||
} else if (lda < MAX(1,m)) {
|
||||
info = 6;
|
||||
} else if (incx == 0) {
|
||||
info = 8;
|
||||
} else if (incy == 0) {
|
||||
info = 11;
|
||||
}
|
||||
if (info != 0) {
|
||||
xerbla('CGEMV ',info);
|
||||
return;
|
||||
}
|
||||
|
||||
// Quick return if possible.
|
||||
|
||||
if ((m == 0) || (n == 0) ||
|
||||
((alpha == zero) && (beta == one))) {
|
||||
return;
|
||||
}
|
||||
|
||||
noconj = lsame(trans,'T');
|
||||
|
||||
// Set LENX and LENY, the lengths of the vectors x and y, and set
|
||||
// up the start points in X and Y.
|
||||
|
||||
if (lsame(trans,'N')) {
|
||||
lenx = n;
|
||||
leny = m;
|
||||
} else {
|
||||
lenx = m;
|
||||
leny = n;
|
||||
}
|
||||
if (incx > 0) {
|
||||
kx = 0;
|
||||
} else {
|
||||
kx = 0 - (lenx-1)*incx;
|
||||
}
|
||||
if (incy > 0) {
|
||||
ky = 0;
|
||||
} else {
|
||||
ky = 0 - (leny-1)*incy;
|
||||
}
|
||||
|
||||
// Start the operations. In this version the elements of A are
|
||||
// accessed sequentially with one pass through A.
|
||||
|
||||
// First form y := beta*y.
|
||||
|
||||
if (beta != one) {
|
||||
if (incy == 1) {
|
||||
if (beta == zero) {
|
||||
for ( i = 0; i < leny; i++) {
|
||||
y[i] = zero;
|
||||
}
|
||||
} else {
|
||||
for ( i = 0; i < leny; i++) {
|
||||
y[i] = beta*y[i];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
iy = ky;
|
||||
if (beta == zero) {
|
||||
for ( i = 0; i < leny; i++) {
|
||||
y[iy] = zero;
|
||||
iy = iy + incy;
|
||||
}
|
||||
} else {
|
||||
for ( i = 0; i < leny; i++) {
|
||||
y[iy] = beta*y[iy];
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if (alpha == zero) {return;}
|
||||
if (lsame(trans,'N')) {
|
||||
|
||||
// Form y := alpha*A*x + y.
|
||||
|
||||
jx = kx;
|
||||
if (incy == 1) {
|
||||
for ( j = 0; j < n; j++) {
|
||||
if (x[jx] != zero) {
|
||||
temp = alpha*x[jx];
|
||||
for ( i = 0; i < m; i++) {
|
||||
y[i] = y[i] + temp*a[i+j*lda];
|
||||
}
|
||||
}
|
||||
jx = jx + incx;
|
||||
}
|
||||
} else {
|
||||
for ( j = 0; j < n; j++) {
|
||||
if (x[jx] != zero) {
|
||||
temp = alpha*x[jx];
|
||||
iy = ky;
|
||||
for ( i = 0; i < m; i++) {
|
||||
y[iy] = y[iy] + temp*a[i+j*lda];
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
jx = jx + incx;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
|
||||
// Form y := alpha*A**T*x + y or y := alpha*A**H*x + y.
|
||||
|
||||
jy = ky;
|
||||
if (incx == 1) {
|
||||
for ( j = 0; j < n; j++) {
|
||||
temp = zero;
|
||||
if (noconj) {
|
||||
for ( i = 0; i < m; i++) {
|
||||
temp = temp + a[i+j*lda]*x[i];
|
||||
}
|
||||
} else {
|
||||
for ( i = 0; i < m; i++) {
|
||||
temp = temp + (~a[i+j*lda])*x[i];
|
||||
}
|
||||
}
|
||||
y[jy] = y[jy] + alpha*temp;
|
||||
jy = jy + incy;
|
||||
}
|
||||
} else {
|
||||
for ( j = 0; j < n; j++) {
|
||||
temp = zero;
|
||||
ix = kx;
|
||||
if (noconj) {
|
||||
for ( i = 0; i < m; i++) {
|
||||
temp = temp + a[i+j*lda]*x[ix];
|
||||
ix = ix + incx;
|
||||
}
|
||||
} else {
|
||||
for ( i = 0; i < m; i++) {
|
||||
temp = temp + (~a[i+j*lda])*x[ix];
|
||||
ix = ix + incx;
|
||||
}
|
||||
}
|
||||
y[jy] = y[jy] + alpha*temp;
|
||||
jy = jy + incy;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of CGEMV .
|
||||
|
||||
}
|
||||
187
src/cgerc.c
187
src/cgerc.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CGERC
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -126,99 +131,89 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CGERC(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
|
||||
*
|
||||
* -- Reference BLAS level2 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
COMPLEX ALPHA
|
||||
INTEGER INCX,INCY,LDA,M,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX A(LDA,*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0E+0,0.0E+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
COMPLEX TEMP
|
||||
INTEGER I,INFO,IX,J,JY,KX
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC CONJG,MAX
|
||||
* ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
IF (M.LT.0) THEN
|
||||
INFO = 1
|
||||
ELSE IF (N.LT.0) THEN
|
||||
INFO = 2
|
||||
ELSE IF (INCX.EQ.0) THEN
|
||||
INFO = 5
|
||||
ELSE IF (INCY.EQ.0) THEN
|
||||
INFO = 7
|
||||
ELSE IF (LDA.LT.MAX(1,M)) THEN
|
||||
INFO = 9
|
||||
END IF
|
||||
IF (INFO.NE.0) THEN
|
||||
CALL XERBLA('CGERC ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF ((M.EQ.0) .OR. (N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
|
||||
*
|
||||
* Start the operations. In this version the elements of A are
|
||||
* accessed sequentially with one pass through A.
|
||||
*
|
||||
IF (INCY.GT.0) THEN
|
||||
JY = 1
|
||||
ELSE
|
||||
JY = 1 - (N-1)*INCY
|
||||
END IF
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 20 J = 1,N
|
||||
IF (Y(JY).NE.ZERO) THEN
|
||||
TEMP = ALPHA*CONJG(Y(JY))
|
||||
DO 10 I = 1,M
|
||||
A(I,J) = A(I,J) + X(I)*TEMP
|
||||
10 CONTINUE
|
||||
END IF
|
||||
JY = JY + INCY
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
IF (INCX.GT.0) THEN
|
||||
KX = 1
|
||||
ELSE
|
||||
KX = 1 - (M-1)*INCX
|
||||
END IF
|
||||
DO 40 J = 1,N
|
||||
IF (Y(JY).NE.ZERO) THEN
|
||||
TEMP = ALPHA*CONJG(Y(JY))
|
||||
IX = KX
|
||||
DO 30 I = 1,M
|
||||
A(I,J) = A(I,J) + X(IX)*TEMP
|
||||
IX = IX + INCX
|
||||
30 CONTINUE
|
||||
END IF
|
||||
JY = JY + INCY
|
||||
40 CONTINUE
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of CGERC
|
||||
*
|
||||
END
|
||||
*/
|
||||
void cgerc(int m, int n, complex float alpha, complex float x[], int incx,
|
||||
complex float y[], int incy, complex float a[], int lda) {
|
||||
|
||||
// -- Reference BLAS level2 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Parameters ..
|
||||
complex float zero = (0.0E+0, 0.0E+0 * I);
|
||||
|
||||
// .. Local Scalars ..
|
||||
complex float temp;
|
||||
int i,info,ix,j,jy,kx;
|
||||
|
||||
// Test the input parameters.
|
||||
|
||||
info = 0;
|
||||
if (m < 0) {
|
||||
info = 1;
|
||||
} else if (n < 0) {
|
||||
info = 2;
|
||||
} else if (incx == 0) {
|
||||
info = 5;
|
||||
} else if (incy == 0) {
|
||||
info = 7;
|
||||
} else if (lda < max(1,m)) {
|
||||
info = 9;
|
||||
}
|
||||
if (info != 0) {
|
||||
xerbla('CGERC ',info);
|
||||
return;
|
||||
}
|
||||
|
||||
// Quick return if possible.
|
||||
|
||||
if ((m == 0) || (n == 0) || (alpha == zero)) return;
|
||||
|
||||
// Start the operations. In this version the elements of A are
|
||||
// accessed sequentially with one pass through A.
|
||||
|
||||
if (incy > 0) {
|
||||
jy = 0;
|
||||
}
|
||||
else {
|
||||
jy = 0 - (n-1)*incy;
|
||||
}
|
||||
if (incx == 1) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if (y[jy] != zero) {
|
||||
temp = alpha*(~y[jy]);
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
a[i+j*lda] = a[i+j*lda] + x[i]*temp;
|
||||
}
|
||||
}
|
||||
jy = jy + incy;
|
||||
}
|
||||
}
|
||||
else {
|
||||
if (incx > 0) {
|
||||
kx = 0;
|
||||
}
|
||||
else {
|
||||
kx = 0 - (m-1)*incx;
|
||||
}
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if (y[jy] != zero) {
|
||||
temp = alpha*(~y[jy]);
|
||||
ix = kx;
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
a[i+j*lda] = a[i+j*lda] + x[ix]*temp;
|
||||
ix = ix + incx;
|
||||
}
|
||||
}
|
||||
jy = jy + incy;
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of CGERC .
|
||||
|
||||
}
|
||||
88
src/cscal.c
88
src/cscal.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CSCAL
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -74,48 +79,41 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CSCAL(N,CA,CX,INCX)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
COMPLEX CA
|
||||
INTEGER INCX,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
INTEGER I,NINCX
|
||||
* ..
|
||||
* .. Parameters ..
|
||||
COMPLEX ONE
|
||||
PARAMETER (ONE= (1.0E+0,0.0E+0))
|
||||
* ..
|
||||
IF (N.LE.0 .OR. INCX.LE.0 .OR. CA.EQ.ONE) RETURN
|
||||
IF (INCX.EQ.1) THEN
|
||||
*
|
||||
* code for increment equal to 1
|
||||
*
|
||||
DO I = 1,N
|
||||
CX(I) = CA*CX(I)
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for increment not equal to 1
|
||||
*
|
||||
NINCX = N*INCX
|
||||
DO I = 1,NINCX,INCX
|
||||
CX(I) = CA*CX(I)
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of CSCAL
|
||||
*
|
||||
END
|
||||
*/
|
||||
void cscal(int n, float complex ca, float complex *cx, int incx) {
|
||||
|
||||
// -- Reference BLAS level1 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
//
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
int i;
|
||||
|
||||
if ( n <= 0 || incx <= 0 ) {
|
||||
return;
|
||||
}
|
||||
|
||||
if (incx == 1) {
|
||||
//
|
||||
// code for increment equal to 1
|
||||
//
|
||||
for (i = 0; i < n; i++) {
|
||||
cx[i] = ca * cx[i];
|
||||
}
|
||||
} else {
|
||||
//
|
||||
// code for increment not equal to 1
|
||||
//
|
||||
for (i = 0; i < n; i++) {
|
||||
cx[i*incx] = ca * cx[i*incx];
|
||||
}
|
||||
}
|
||||
return;
|
||||
|
||||
// End of CSCAL
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
120
src/csrot.c
120
src/csrot.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CSROT
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -94,60 +99,61 @@
|
|||
*> \ingroup complex_blas_level1
|
||||
*
|
||||
* =====================================================================
|
||||
SUBROUTINE CSROT( N, CX, INCX, CY, INCY, C, S )
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX, INCY, N
|
||||
REAL C, S
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX( * ), CY( * )
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
INTEGER I, IX, IY
|
||||
COMPLEX CTEMP
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
IF( N.LE.0 )
|
||||
$ RETURN
|
||||
IF( INCX.EQ.1 .AND. INCY.EQ.1 ) THEN
|
||||
*
|
||||
* code for both increments equal to 1
|
||||
*
|
||||
DO I = 1, N
|
||||
CTEMP = C*CX( I ) + S*CY( I )
|
||||
CY( I ) = C*CY( I ) - S*CX( I )
|
||||
CX( I ) = CTEMP
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for unequal increments or equal increments not equal
|
||||
* to 1
|
||||
*
|
||||
IX = 1
|
||||
IY = 1
|
||||
IF( INCX.LT.0 )
|
||||
$ IX = ( -N+1 )*INCX + 1
|
||||
IF( INCY.LT.0 )
|
||||
$ IY = ( -N+1 )*INCY + 1
|
||||
DO I = 1, N
|
||||
CTEMP = C*CX( IX ) + S*CY( IY )
|
||||
CY( IY ) = C*CY( IY ) - S*CX( IX )
|
||||
CX( IX ) = CTEMP
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of CSROT
|
||||
*
|
||||
END
|
||||
*/
|
||||
void csrot ( int n, float complex cx[], int incx, float complex cy[],
|
||||
int incy, float c, float s ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
float complex ctemp;
|
||||
int i, ix, iy;
|
||||
|
||||
if ( n <= 0 ) {
|
||||
return;
|
||||
}
|
||||
|
||||
if ( incx == 1 && incy == 1 ) {
|
||||
|
||||
// code for both increments equal to 1
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
ctemp = c * cx[i] + s * cy[i];
|
||||
cy[i] = c * cy[i] - s * cx[i];
|
||||
cx[i] = ctemp;
|
||||
}
|
||||
} else {
|
||||
|
||||
// code for unequal increments or equal increments not equal
|
||||
// to 1
|
||||
|
||||
if ( 0 <= incx ) {
|
||||
ix = 0;
|
||||
} else {
|
||||
ix = ( -n + 1 ) * incx;
|
||||
}
|
||||
|
||||
if ( 0 <= incy ) {
|
||||
iy = 0;
|
||||
} else {
|
||||
iy = ( -n + 1 ) * incy;
|
||||
}
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
ctemp = c * cx[ix] + s * cy[iy];
|
||||
cy[iy] = c * cy[iy] - s * cx[ix];
|
||||
cx[ix] = ctemp;
|
||||
ix = ix + incx;
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of CSROT
|
||||
|
||||
}
|
||||
|
|
|
|||
90
src/csscal.c
90
src/csscal.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CSSCAL
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -74,51 +79,40 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CSSCAL(N,SA,CX,INCX)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
REAL SA
|
||||
INTEGER INCX,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
INTEGER I,NINCX
|
||||
* ..
|
||||
* .. Parameters ..
|
||||
REAL ONE
|
||||
PARAMETER (ONE=1.0E+0)
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC AIMAG,CMPLX,REAL
|
||||
* ..
|
||||
IF (N.LE.0 .OR. INCX.LE.0 .OR. SA.EQ.ONE) RETURN
|
||||
IF (INCX.EQ.1) THEN
|
||||
*
|
||||
* code for increment equal to 1
|
||||
*
|
||||
DO I = 1,N
|
||||
CX(I) = CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I)))
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for increment not equal to 1
|
||||
*
|
||||
NINCX = N*INCX
|
||||
DO I = 1,NINCX,INCX
|
||||
CX(I) = CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I)))
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of CSSCAL
|
||||
*
|
||||
END
|
||||
*/
|
||||
void csscal ( int n, float sa, float complex cx[], int incx ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
int i;
|
||||
|
||||
if ( n <= 0 || incx <= 0 ) {
|
||||
return;
|
||||
}
|
||||
|
||||
if ( incx == 1 ) {
|
||||
|
||||
// code for increment equal to 1
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
cx[i] = sa * cx[i];
|
||||
}
|
||||
} else {
|
||||
|
||||
// code for increment not equal to 1
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
cx[i*incx] = sa * cx[i*incx];
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of CSSCAL
|
||||
|
||||
}
|
||||
|
|
|
|||
108
src/cswap.c
108
src/cswap.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CSWAP
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -77,53 +82,56 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CSWAP(N,CX,INCX,CY,INCY)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*),CY(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
COMPLEX CTEMP
|
||||
INTEGER I,IX,IY
|
||||
* ..
|
||||
IF (N.LE.0) RETURN
|
||||
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
|
||||
*
|
||||
* code for both increments equal to 1
|
||||
DO I = 1,N
|
||||
CTEMP = CX(I)
|
||||
CX(I) = CY(I)
|
||||
CY(I) = CTEMP
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for unequal increments or equal increments not equal
|
||||
* to 1
|
||||
*
|
||||
IX = 1
|
||||
IY = 1
|
||||
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
|
||||
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
|
||||
DO I = 1,N
|
||||
CTEMP = CX(IX)
|
||||
CX(IX) = CY(IY)
|
||||
CY(IY) = CTEMP
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of CSWAP
|
||||
*
|
||||
END
|
||||
*/
|
||||
void cswap ( int n, float complex cx[], int incx, float complex cy[],
|
||||
int incy ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
float complex ctemp;
|
||||
int i, ix, iy;
|
||||
|
||||
if ( n <= 0 ) return;
|
||||
|
||||
if ( incx == 1 && incy == 1 ) {
|
||||
// code for both increments equal to 1
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
ctemp = cx[i];
|
||||
cx[i] = cy[i];
|
||||
cy[i] = ctemp;
|
||||
}
|
||||
} else {
|
||||
|
||||
// code for unequal increments or equal increments not equal
|
||||
// to 1
|
||||
|
||||
if ( 0 <= incx ) {
|
||||
ix = 0;
|
||||
} else {
|
||||
ix = ( -n + 1 ) * incx;
|
||||
}
|
||||
|
||||
if ( 0 <= incy ) {
|
||||
iy = 0;
|
||||
} else {
|
||||
iy = ( -n + 1 ) * incy;
|
||||
}
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
ctemp = cx[ix];
|
||||
cx[ix] = cy[iy];
|
||||
cy[iy] = ctemp;
|
||||
ix = ix + incx;
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of CSWAP
|
||||
|
||||
}
|
||||
|
|
|
|||
459
src/ctrmv.c
459
src/ctrmv.c
|
|
@ -1,3 +1,9 @@
|
|||
#include <complex.h>
|
||||
#include <stdbool.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b CTRMV
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -143,228 +149,231 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE CTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
|
||||
*
|
||||
* -- Reference BLAS level2 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,LDA,N
|
||||
CHARACTER DIAG,TRANS,UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX A(LDA,*),X(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0E+0,0.0E+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
COMPLEX TEMP
|
||||
INTEGER I,INFO,IX,J,JX,KX
|
||||
LOGICAL NOCONJ,NOUNIT
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC CONJG,MAX
|
||||
* ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
|
||||
INFO = 1
|
||||
ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
|
||||
+ .NOT.LSAME(TRANS,'C')) THEN
|
||||
INFO = 2
|
||||
ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
|
||||
INFO = 3
|
||||
ELSE IF (N.LT.0) THEN
|
||||
INFO = 4
|
||||
ELSE IF (LDA.LT.MAX(1,N)) THEN
|
||||
INFO = 6
|
||||
ELSE IF (INCX.EQ.0) THEN
|
||||
INFO = 8
|
||||
END IF
|
||||
IF (INFO.NE.0) THEN
|
||||
CALL XERBLA('CTRMV ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF (N.EQ.0) RETURN
|
||||
*
|
||||
NOCONJ = LSAME(TRANS,'T')
|
||||
NOUNIT = LSAME(DIAG,'N')
|
||||
*
|
||||
* Set up the start point in X if the increment is not unity. This
|
||||
* will be ( N - 1 )*INCX too small for descending loops.
|
||||
*
|
||||
IF (INCX.LE.0) THEN
|
||||
KX = 1 - (N-1)*INCX
|
||||
ELSE IF (INCX.NE.1) THEN
|
||||
KX = 1
|
||||
END IF
|
||||
*
|
||||
* Start the operations. In this version the elements of A are
|
||||
* accessed sequentially with one pass through A.
|
||||
*
|
||||
IF (LSAME(TRANS,'N')) THEN
|
||||
*
|
||||
* Form x := A*x.
|
||||
*
|
||||
IF (LSAME(UPLO,'U')) THEN
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 20 J = 1,N
|
||||
IF (X(J).NE.ZERO) THEN
|
||||
TEMP = X(J)
|
||||
DO 10 I = 1,J - 1
|
||||
X(I) = X(I) + TEMP*A(I,J)
|
||||
10 CONTINUE
|
||||
IF (NOUNIT) X(J) = X(J)*A(J,J)
|
||||
END IF
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 40 J = 1,N
|
||||
IF (X(JX).NE.ZERO) THEN
|
||||
TEMP = X(JX)
|
||||
IX = KX
|
||||
DO 30 I = 1,J - 1
|
||||
X(IX) = X(IX) + TEMP*A(I,J)
|
||||
IX = IX + INCX
|
||||
30 CONTINUE
|
||||
IF (NOUNIT) X(JX) = X(JX)*A(J,J)
|
||||
END IF
|
||||
JX = JX + INCX
|
||||
40 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 60 J = N,1,-1
|
||||
IF (X(J).NE.ZERO) THEN
|
||||
TEMP = X(J)
|
||||
DO 50 I = N,J + 1,-1
|
||||
X(I) = X(I) + TEMP*A(I,J)
|
||||
50 CONTINUE
|
||||
IF (NOUNIT) X(J) = X(J)*A(J,J)
|
||||
END IF
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
KX = KX + (N-1)*INCX
|
||||
JX = KX
|
||||
DO 80 J = N,1,-1
|
||||
IF (X(JX).NE.ZERO) THEN
|
||||
TEMP = X(JX)
|
||||
IX = KX
|
||||
DO 70 I = N,J + 1,-1
|
||||
X(IX) = X(IX) + TEMP*A(I,J)
|
||||
IX = IX - INCX
|
||||
70 CONTINUE
|
||||
IF (NOUNIT) X(JX) = X(JX)*A(J,J)
|
||||
END IF
|
||||
JX = JX - INCX
|
||||
80 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
ELSE
|
||||
*
|
||||
* Form x := A**T*x or x := A**H*x.
|
||||
*
|
||||
IF (LSAME(UPLO,'U')) THEN
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 110 J = N,1,-1
|
||||
TEMP = X(J)
|
||||
IF (NOCONJ) THEN
|
||||
IF (NOUNIT) TEMP = TEMP*A(J,J)
|
||||
DO 90 I = J - 1,1,-1
|
||||
TEMP = TEMP + A(I,J)*X(I)
|
||||
90 CONTINUE
|
||||
ELSE
|
||||
IF (NOUNIT) TEMP = TEMP*CONJG(A(J,J))
|
||||
DO 100 I = J - 1,1,-1
|
||||
TEMP = TEMP + CONJG(A(I,J))*X(I)
|
||||
100 CONTINUE
|
||||
END IF
|
||||
X(J) = TEMP
|
||||
110 CONTINUE
|
||||
ELSE
|
||||
JX = KX + (N-1)*INCX
|
||||
DO 140 J = N,1,-1
|
||||
TEMP = X(JX)
|
||||
IX = JX
|
||||
IF (NOCONJ) THEN
|
||||
IF (NOUNIT) TEMP = TEMP*A(J,J)
|
||||
DO 120 I = J - 1,1,-1
|
||||
IX = IX - INCX
|
||||
TEMP = TEMP + A(I,J)*X(IX)
|
||||
120 CONTINUE
|
||||
ELSE
|
||||
IF (NOUNIT) TEMP = TEMP*CONJG(A(J,J))
|
||||
DO 130 I = J - 1,1,-1
|
||||
IX = IX - INCX
|
||||
TEMP = TEMP + CONJG(A(I,J))*X(IX)
|
||||
130 CONTINUE
|
||||
END IF
|
||||
X(JX) = TEMP
|
||||
JX = JX - INCX
|
||||
140 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 170 J = 1,N
|
||||
TEMP = X(J)
|
||||
IF (NOCONJ) THEN
|
||||
IF (NOUNIT) TEMP = TEMP*A(J,J)
|
||||
DO 150 I = J + 1,N
|
||||
TEMP = TEMP + A(I,J)*X(I)
|
||||
150 CONTINUE
|
||||
ELSE
|
||||
IF (NOUNIT) TEMP = TEMP*CONJG(A(J,J))
|
||||
DO 160 I = J + 1,N
|
||||
TEMP = TEMP + CONJG(A(I,J))*X(I)
|
||||
160 CONTINUE
|
||||
END IF
|
||||
X(J) = TEMP
|
||||
170 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 200 J = 1,N
|
||||
TEMP = X(JX)
|
||||
IX = JX
|
||||
IF (NOCONJ) THEN
|
||||
IF (NOUNIT) TEMP = TEMP*A(J,J)
|
||||
DO 180 I = J + 1,N
|
||||
IX = IX + INCX
|
||||
TEMP = TEMP + A(I,J)*X(IX)
|
||||
180 CONTINUE
|
||||
ELSE
|
||||
IF (NOUNIT) TEMP = TEMP*CONJG(A(J,J))
|
||||
DO 190 I = J + 1,N
|
||||
IX = IX + INCX
|
||||
TEMP = TEMP + CONJG(A(I,J))*X(IX)
|
||||
190 CONTINUE
|
||||
END IF
|
||||
X(JX) = TEMP
|
||||
JX = JX + INCX
|
||||
200 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of CTRMV
|
||||
*
|
||||
END
|
||||
*/
|
||||
void ctrmv(char uplo, char trans, char diag, int n, complex float a[], int lda, complex float x[], int incx) {
|
||||
|
||||
// -- Reference BLAS level2 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Parameters ..
|
||||
complex float zero = (0.0E+0,0.0E+0*I);
|
||||
// ..
|
||||
// .. Local Scalars ..
|
||||
complex float temp;
|
||||
int i,info,ix,j,jx,kx;
|
||||
bool noconj,nounit;
|
||||
|
||||
// Test the input parameters.
|
||||
//
|
||||
info = 0;
|
||||
if ( !lsame( uplo, 'U' ) && !lsame( uplo, 'L' ) ) {
|
||||
info = 1;
|
||||
} else if ( !lsame( trans, 'N' ) && !lsame( trans, 'T' ) &&
|
||||
!lsame( trans, 'C' ) ) {
|
||||
info = 2;
|
||||
} else if ( !lsame( diag, 'U' ) && !lsame( diag, 'N' ) ) {
|
||||
info = 3;
|
||||
} else if (n < 0) {
|
||||
info = 4;
|
||||
} else if ( lda < i4_max( 1, n ) ) {
|
||||
info = 6;
|
||||
} else if ( incx == 0 ) {
|
||||
info = 8;
|
||||
}
|
||||
if ( info != 0 ) {
|
||||
xerbla('CTRMV', info);
|
||||
return;
|
||||
}
|
||||
//
|
||||
// Quick return if possible.
|
||||
//
|
||||
if ( n == 0 ) {
|
||||
return;
|
||||
}
|
||||
//
|
||||
noconj = lsame( trans, 'T' );
|
||||
nounit = lsame( diag, 'N' );
|
||||
|
||||
// Set up the start point in X if the increment is not unity. This
|
||||
// will be ( N - 1 )*INCX too small for descending loops.
|
||||
|
||||
if ( incx <= 0 ) {
|
||||
kx = 0 - ( n-1 ) * incx;
|
||||
} else if ( incx != 1 ) {
|
||||
kx = 0;
|
||||
}
|
||||
|
||||
// Start the operations. In this version the elements of A are
|
||||
// accessed sequentially with one pass through A.
|
||||
|
||||
if ( lsame( trans, 'N' ) ) {
|
||||
|
||||
// Form x := A*x.
|
||||
|
||||
if ( lsame(uplo,'U') ) {
|
||||
if (incx == 1) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( x[j] != zero ) {
|
||||
temp = x[j];
|
||||
for ( i = 0; i < j - 1; i++ ) {
|
||||
x[i] = x[i] + temp * a[i+j*lda];
|
||||
}
|
||||
if ( nounit ) {
|
||||
x[j] = x[j] * a[j+j*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
jx = kx;
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if (x[jx] != zero) {
|
||||
temp = x[jx];
|
||||
ix = kx;
|
||||
for ( i = 0; i < j; i++ ) {
|
||||
x[ix] = x[ix] + temp * a[i+j*lda];
|
||||
ix = ix + incx;
|
||||
}
|
||||
if ( nounit ) {
|
||||
x[jx] = x[jx] * a[j+j*lda];
|
||||
}
|
||||
}
|
||||
jx = jx + incx;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if ( incx == 1 ) {
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
if ( x[j] != zero) {
|
||||
temp = x[j];
|
||||
for ( i = n - 1; j < i; i-- ) {
|
||||
x[i] = x[i] + temp * a[i+j*lda];
|
||||
}
|
||||
if ( nounit ) {
|
||||
x[j] = x[j] * a[j+j*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
kx = kx + ( n - 1 ) * incx;
|
||||
jx = kx;
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
if (x[jx] != zero) {
|
||||
temp = x[jx];
|
||||
ix = kx;
|
||||
for ( i = n - 1; j < i; i-- ) {
|
||||
x[ix] = x[ix] + temp * a[i+j*lda];
|
||||
ix = ix - incx;
|
||||
}
|
||||
if ( nounit ) {
|
||||
x[jx] = x[jx] * a[j+j*lda];
|
||||
}
|
||||
}
|
||||
jx = jx - incx;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
|
||||
// Form x := A**T*x or x := A**H*x.
|
||||
|
||||
if ( lsame(uplo,'U') ) {
|
||||
if ( incx == 1 ) {
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
temp = x[j];
|
||||
if ( noconj ) {
|
||||
if ( nounit ) {
|
||||
temp = temp * a[j+j*lda];
|
||||
}
|
||||
for ( i = j - 1; 0 <= i; i-- ) {
|
||||
temp = temp + a[i+j*lda] * x[i];
|
||||
}
|
||||
} else {
|
||||
if (nounit) {
|
||||
temp = temp * (~a[j+j*lda]);
|
||||
}
|
||||
for ( i = j-1; 0 <= i; i-- ) {
|
||||
temp = temp + (~a[i+j*lda])*x[i];
|
||||
}
|
||||
}
|
||||
x[j] = temp;
|
||||
}
|
||||
} else {
|
||||
jx = kx + ( n-1 ) * incx;
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
temp = x[jx];
|
||||
ix = jx;
|
||||
if ( noconj ) {
|
||||
if ( nounit ) {
|
||||
temp = temp * a[j+j*lda];
|
||||
}
|
||||
for ( i = j - 1; 0 <= i; i-- ) {
|
||||
ix = ix - incx;
|
||||
temp = temp + a[i+j*lda] * x[ix];
|
||||
}
|
||||
} else {
|
||||
if (nounit) {
|
||||
temp = temp * (~a[j+j*lda]);
|
||||
}
|
||||
for ( i = j-1; 0 <= i; i-- ) {
|
||||
ix = ix - incx;
|
||||
temp = temp + (~a[i+j*lda]) * x[ix];
|
||||
}
|
||||
}
|
||||
x[jx] = temp;
|
||||
jx = jx - incx;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if ( incx == 1 ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
temp = x[j];
|
||||
if ( noconj ) {
|
||||
if (nounit) {
|
||||
temp = temp*a[j+j*lda];
|
||||
}
|
||||
for ( i = j + 1; i < n; i++ ) {
|
||||
temp = temp + a[i+j*lda] * x[i];
|
||||
}
|
||||
} else {
|
||||
if (nounit) temp = temp * (~a[j+j*lda]);
|
||||
for ( i = j + 1; i < n; i++ ) {
|
||||
temp = temp + (~a[i+j*lda]) * x[i];
|
||||
}
|
||||
}
|
||||
x[j] = temp;
|
||||
}
|
||||
} else {
|
||||
jx = kx;
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
temp = x[jx];
|
||||
ix = jx;
|
||||
if ( noconj ) {
|
||||
if ( nounit ) {
|
||||
temp = temp*a[j+j*lda];
|
||||
}
|
||||
for ( i = j + 1; i < n; i++ ) {
|
||||
ix = ix + incx;
|
||||
temp = temp + a[i+j*lda] * x[ix];
|
||||
}
|
||||
} else {
|
||||
if (nounit) temp = temp * (~a[j+j*lda]);
|
||||
for ( i = j + 1; i < n; i++ ) {
|
||||
ix = ix + incx;
|
||||
temp = temp + (~a[i+j*lda]) * x[ix];
|
||||
}
|
||||
}
|
||||
x[jx] = temp;
|
||||
jx = jx + incx;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of CTRMV .
|
||||
|
||||
}
|
||||
|
|
|
|||
129
src/ddot.c
129
src/ddot.c
|
|
@ -1,3 +1,4 @@
|
|||
/*
|
||||
*> \brief \b DDOT
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -78,71 +79,63 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE PRECISION DX(*),DY(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
DOUBLE PRECISION DTEMP
|
||||
INTEGER I,IX,IY,M,MP1
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MOD
|
||||
* ..
|
||||
DDOT = 0.0d0
|
||||
DTEMP = 0.0d0
|
||||
IF (N.LE.0) RETURN
|
||||
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
|
||||
*
|
||||
* code for both increments equal to 1
|
||||
*
|
||||
*
|
||||
* clean-up loop
|
||||
*
|
||||
M = MOD(N,5)
|
||||
IF (M.NE.0) THEN
|
||||
DO I = 1,M
|
||||
DTEMP = DTEMP + DX(I)*DY(I)
|
||||
END DO
|
||||
IF (N.LT.5) THEN
|
||||
DDOT=DTEMP
|
||||
RETURN
|
||||
END IF
|
||||
END IF
|
||||
MP1 = M + 1
|
||||
DO I = MP1,N,5
|
||||
DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) +
|
||||
$ DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for unequal increments or equal increments
|
||||
* not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
IY = 1
|
||||
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
|
||||
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
|
||||
DO I = 1,N
|
||||
DTEMP = DTEMP + DX(IX)*DY(IY)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
END DO
|
||||
END IF
|
||||
DDOT = DTEMP
|
||||
RETURN
|
||||
*
|
||||
* End of DDOT
|
||||
*
|
||||
END
|
||||
*/
|
||||
double ddot(int n, double dx[], int incx, double dy[], int incy) {
|
||||
|
||||
// -- Reference BLAS level1 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
|
||||
// =====================================================================
|
||||
//
|
||||
// .. Local Scalars ..
|
||||
double dtemp;
|
||||
int i, ix, iy, m;
|
||||
|
||||
// ..
|
||||
dtemp = 0;
|
||||
if (n <= 0) return dtemp;
|
||||
if (incx == 1 && incy == 1) {
|
||||
//
|
||||
// code for both increments equal to 1
|
||||
//
|
||||
//
|
||||
// clean-up loop
|
||||
//
|
||||
m = n % 5;
|
||||
if (m != 0) {
|
||||
for (i = 0; i < m; i++) {
|
||||
dtemp = dtemp + dx[i]*dy[i];
|
||||
}
|
||||
if (n < 5) return dtemp;
|
||||
}
|
||||
// mp1 = m + 1;
|
||||
for (i = m; i < n; i += 5) {
|
||||
dtemp += dx[i ]*dy[i ] +
|
||||
dx[i+1]*dy[i+1] +
|
||||
dx[i+2]*dy[i+2] +
|
||||
dx[i+3]*dy[i+3] +
|
||||
dx[i+4]*dy[i+4];
|
||||
}
|
||||
} else {
|
||||
//
|
||||
// code for unequal increments or equal increments
|
||||
// not equal to 1
|
||||
//
|
||||
ix = 0;
|
||||
iy = 0;
|
||||
if (incx < 0) ix = (-n+1)*incx;
|
||||
if (incy < 0) iy = (-n+1)*incy;
|
||||
for (i = 0; i < n; i++) {
|
||||
dtemp += dx[ix]*dy[iy];
|
||||
ix += incx;
|
||||
iy += incy;
|
||||
}
|
||||
}
|
||||
|
||||
return dtemp;
|
||||
|
||||
// End DDOT
|
||||
|
||||
}
|
||||
|
|
|
|||
120
src/icamax.c
120
src/icamax.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b ICAMAX
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -67,61 +72,60 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
INTEGER FUNCTION ICAMAX(N,CX,INCX)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
REAL SMAX
|
||||
INTEGER I,IX
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
REAL SCABS1
|
||||
EXTERNAL SCABS1
|
||||
* ..
|
||||
ICAMAX = 0
|
||||
IF (N.LT.1 .OR. INCX.LE.0) RETURN
|
||||
ICAMAX = 1
|
||||
IF (N.EQ.1) RETURN
|
||||
IF (INCX.EQ.1) THEN
|
||||
*
|
||||
* code for increment equal to 1
|
||||
*
|
||||
SMAX = SCABS1(CX(1))
|
||||
DO I = 2,N
|
||||
IF (SCABS1(CX(I)).GT.SMAX) THEN
|
||||
ICAMAX = I
|
||||
SMAX = SCABS1(CX(I))
|
||||
END IF
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for increment not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
SMAX = SCABS1(CX(1))
|
||||
IX = IX + INCX
|
||||
DO I = 2,N
|
||||
IF (SCABS1(CX(IX)).GT.SMAX) THEN
|
||||
ICAMAX = I
|
||||
SMAX = SCABS1(CX(IX))
|
||||
END IF
|
||||
IX = IX + INCX
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of ICAMAX
|
||||
*
|
||||
END
|
||||
*/
|
||||
int icamax ( int n, float complex x[], int incx ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
float smax;
|
||||
int i, ix, value;
|
||||
|
||||
value = 0;
|
||||
|
||||
if ( n < 1 || incx <= 0 ) {
|
||||
return value;
|
||||
}
|
||||
|
||||
value = 1;
|
||||
|
||||
if ( n == 1 ) {
|
||||
return value;
|
||||
}
|
||||
|
||||
if ( incx != 1 ) {
|
||||
ix = 0;
|
||||
smax = cabs1 ( x[0] );
|
||||
ix = ix + incx;
|
||||
|
||||
for ( i = 1; i < n; i++ )
|
||||
{
|
||||
if ( smax < cabs1 ( x[ix] ) )
|
||||
{
|
||||
value = i + 1;
|
||||
smax = cabs1 ( x[ix] );
|
||||
}
|
||||
ix = ix + incx;
|
||||
}
|
||||
}
|
||||
else {
|
||||
smax = cabs1 ( x[0] );
|
||||
for ( i = 1; i < n; i++ )
|
||||
{
|
||||
if ( smax < cabs1 ( x[i] ) )
|
||||
{
|
||||
value = i + 1;
|
||||
smax = cabs1 ( x[i] );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return value;
|
||||
|
||||
// End of ICAMAX
|
||||
|
||||
}
|
||||
|
|
|
|||
115
src/isamax.c
115
src/isamax.c
|
|
@ -1,3 +1,6 @@
|
|||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b ISAMAX
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -67,60 +70,58 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
INTEGER FUNCTION ISAMAX(N,SX,INCX)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL SX(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
REAL SMAX
|
||||
INTEGER I,IX
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS
|
||||
* ..
|
||||
ISAMAX = 0
|
||||
IF (N.LT.1 .OR. INCX.LE.0) RETURN
|
||||
ISAMAX = 1
|
||||
IF (N.EQ.1) RETURN
|
||||
IF (INCX.EQ.1) THEN
|
||||
*
|
||||
* code for increment equal to 1
|
||||
*
|
||||
SMAX = ABS(SX(1))
|
||||
DO I = 2,N
|
||||
IF (ABS(SX(I)).GT.SMAX) THEN
|
||||
ISAMAX = I
|
||||
SMAX = ABS(SX(I))
|
||||
END IF
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for increment not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
SMAX = ABS(SX(1))
|
||||
IX = IX + INCX
|
||||
DO I = 2,N
|
||||
IF (ABS(SX(IX)).GT.SMAX) THEN
|
||||
ISAMAX = I
|
||||
SMAX = ABS(SX(IX))
|
||||
END IF
|
||||
IX = IX + INCX
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of ISAMAX
|
||||
*
|
||||
END
|
||||
*/
|
||||
int isamax ( int n, float dx[], int incx ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
float dmax;
|
||||
int i, ix, value;
|
||||
|
||||
value = 0;
|
||||
|
||||
if ( n < 1 || incx <= 0 ) {
|
||||
return value;
|
||||
}
|
||||
|
||||
value = 1;
|
||||
|
||||
if ( n == 1 ) {
|
||||
return value;
|
||||
}
|
||||
|
||||
if ( incx == 1 ) {
|
||||
dmax =fabs ( dx[0] );
|
||||
|
||||
for ( i = 1; i < n; i++ ) {
|
||||
if ( dmax < fabs ( dx[i] ) ) {
|
||||
value = i + 1;
|
||||
dmax = fabs ( dx[i] );
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
ix = 0;
|
||||
dmax = fabs ( dx[0] );
|
||||
ix = ix + incx;
|
||||
|
||||
for ( i = 1; i < n; i++ ) {
|
||||
if ( dmax < fabs ( dx[ix] ) ) {
|
||||
value = i + 1;
|
||||
dmax = fabs ( dx[ix] );
|
||||
}
|
||||
ix = ix + incx;
|
||||
}
|
||||
}
|
||||
|
||||
return value;
|
||||
|
||||
// End of ISAMAX
|
||||
|
||||
}
|
||||
110
src/izamax.c
110
src/izamax.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b IZAMAX
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -67,61 +72,50 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
INTEGER FUNCTION IZAMAX(N,ZX,INCX)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX*16 ZX(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
DOUBLE PRECISION DMAX
|
||||
INTEGER I,IX
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
DOUBLE PRECISION DCABS1
|
||||
EXTERNAL DCABS1
|
||||
* ..
|
||||
IZAMAX = 0
|
||||
IF (N.LT.1 .OR. INCX.LE.0) RETURN
|
||||
IZAMAX = 1
|
||||
IF (N.EQ.1) RETURN
|
||||
IF (INCX.EQ.1) THEN
|
||||
*
|
||||
* code for increment equal to 1
|
||||
*
|
||||
DMAX = DCABS1(ZX(1))
|
||||
DO I = 2,N
|
||||
IF (DCABS1(ZX(I)).GT.DMAX) THEN
|
||||
IZAMAX = I
|
||||
DMAX = DCABS1(ZX(I))
|
||||
END IF
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for increment not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
DMAX = DCABS1(ZX(1))
|
||||
IX = IX + INCX
|
||||
DO I = 2,N
|
||||
IF (DCABS1(ZX(IX)).GT.DMAX) THEN
|
||||
IZAMAX = I
|
||||
DMAX = DCABS1(ZX(IX))
|
||||
END IF
|
||||
IX = IX + INCX
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of IZAMAX
|
||||
*
|
||||
END
|
||||
*/
|
||||
int izamax(int n, complex double zx[], int incx) {
|
||||
|
||||
// -- Reference BLAS level1 routine (version 3.4.0) --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
// November 2011
|
||||
|
||||
// .. Local Scalars ..
|
||||
double dmax;
|
||||
int i, ix, ret_val;
|
||||
|
||||
ret_val = 0;
|
||||
if (n < 1 || incx <= 0) return ret_val;
|
||||
ret_val = 1;
|
||||
if (n == 1) return ret_val;
|
||||
if (incx == 1) {
|
||||
|
||||
// code for increment equal to 1
|
||||
dmax = dcabs1(zx[1]);
|
||||
for (i = 1; i < n; i++) {
|
||||
if (dcabs1(zx[i]) > dmax) {
|
||||
ret_val = i;
|
||||
dmax = dcabs1(zx[i]);
|
||||
}
|
||||
}
|
||||
} else {
|
||||
|
||||
// code for increment not equal to 1
|
||||
|
||||
ix = 1;
|
||||
dmax = dcabs1(zx[1]);
|
||||
ix += incx;
|
||||
for (i = 1; i < n; i++) {
|
||||
if (dcabs1(zx[ix]) > dmax) {
|
||||
ret_val = i;
|
||||
dmax = dcabs1(zx[ix]);
|
||||
}
|
||||
ix += incx;
|
||||
}
|
||||
}
|
||||
|
||||
return ret_val;
|
||||
|
||||
// End of IZAMAX
|
||||
|
||||
}
|
||||
|
|
|
|||
104
src/lsame.c
104
src/lsame.c
|
|
@ -1,3 +1,6 @@
|
|||
#include <stdbool.h>
|
||||
|
||||
/*
|
||||
*> \brief \b LSAME
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -49,74 +52,33 @@
|
|||
*> \ingroup aux_blas
|
||||
*
|
||||
* =====================================================================
|
||||
LOGICAL FUNCTION LSAME(CA,CB)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
CHARACTER CA,CB
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ICHAR
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
INTEGER INTA,INTB,ZCODE
|
||||
* ..
|
||||
*
|
||||
* Test if the characters are equal
|
||||
*
|
||||
LSAME = CA .EQ. CB
|
||||
IF (LSAME) RETURN
|
||||
*
|
||||
* Now test for equivalence if both characters are alphabetic.
|
||||
*
|
||||
ZCODE = ICHAR('Z')
|
||||
*
|
||||
* Use 'Z' rather than 'A' so that ASCII can be detected on Prime
|
||||
* machines, on which ICHAR returns a value with bit 8 set.
|
||||
* ICHAR('A') on Prime machines returns 193 which is the same as
|
||||
* ICHAR('A') on an EBCDIC machine.
|
||||
*
|
||||
INTA = ICHAR(CA)
|
||||
INTB = ICHAR(CB)
|
||||
*
|
||||
IF (ZCODE.EQ.90 .OR. ZCODE.EQ.122) THEN
|
||||
*
|
||||
* ASCII is assumed - ZCODE is the ASCII code of either lower or
|
||||
* upper case 'Z'.
|
||||
*
|
||||
IF (INTA.GE.97 .AND. INTA.LE.122) INTA = INTA - 32
|
||||
IF (INTB.GE.97 .AND. INTB.LE.122) INTB = INTB - 32
|
||||
*
|
||||
ELSE IF (ZCODE.EQ.233 .OR. ZCODE.EQ.169) THEN
|
||||
*
|
||||
* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
|
||||
* upper case 'Z'.
|
||||
*
|
||||
IF (INTA.GE.129 .AND. INTA.LE.137 .OR.
|
||||
+ INTA.GE.145 .AND. INTA.LE.153 .OR.
|
||||
+ INTA.GE.162 .AND. INTA.LE.169) INTA = INTA + 64
|
||||
IF (INTB.GE.129 .AND. INTB.LE.137 .OR.
|
||||
+ INTB.GE.145 .AND. INTB.LE.153 .OR.
|
||||
+ INTB.GE.162 .AND. INTB.LE.169) INTB = INTB + 64
|
||||
*
|
||||
ELSE IF (ZCODE.EQ.218 .OR. ZCODE.EQ.250) THEN
|
||||
*
|
||||
* ASCII is assumed, on Prime machines - ZCODE is the ASCII code
|
||||
* plus 128 of either lower or upper case 'Z'.
|
||||
*
|
||||
IF (INTA.GE.225 .AND. INTA.LE.250) INTA = INTA - 32
|
||||
IF (INTB.GE.225 .AND. INTB.LE.250) INTB = INTB - 32
|
||||
END IF
|
||||
LSAME = INTA .EQ. INTB
|
||||
*
|
||||
* RETURN
|
||||
*
|
||||
* End of LSAME
|
||||
*
|
||||
END
|
||||
*/
|
||||
bool lsame ( char ca, char cb ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// Test if the characters are equal
|
||||
|
||||
if ( ca == cb ) {
|
||||
return true;
|
||||
}
|
||||
|
||||
if ( 'A' <= ca && ca <= 'Z' ) {
|
||||
if ( ca - 'A' == cb - 'a' ) {
|
||||
return true;
|
||||
}
|
||||
} else if ( 'a' <= ca && ca <= 'z' ) {
|
||||
if ( ca - 'a' == cb - 'A' ) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
return false;
|
||||
|
||||
// End of LSAME
|
||||
|
||||
}
|
||||
88
src/sasum.c
88
src/sasum.c
|
|
@ -1,3 +1,6 @@
|
|||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b SASUM
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -68,65 +71,26 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
REAL FUNCTION SASUM(N,SX,INCX)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL SX(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
REAL STEMP
|
||||
INTEGER I,M,MP1,NINCX
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS,MOD
|
||||
* ..
|
||||
SASUM = 0.0e0
|
||||
STEMP = 0.0e0
|
||||
IF (N.LE.0 .OR. INCX.LE.0) RETURN
|
||||
IF (INCX.EQ.1) THEN
|
||||
* code for increment equal to 1
|
||||
*
|
||||
*
|
||||
* clean-up loop
|
||||
*
|
||||
M = MOD(N,6)
|
||||
IF (M.NE.0) THEN
|
||||
DO I = 1,M
|
||||
STEMP = STEMP + ABS(SX(I))
|
||||
END DO
|
||||
IF (N.LT.6) THEN
|
||||
SASUM = STEMP
|
||||
RETURN
|
||||
END IF
|
||||
END IF
|
||||
MP1 = M + 1
|
||||
DO I = MP1,N,6
|
||||
STEMP = STEMP + ABS(SX(I)) + ABS(SX(I+1)) +
|
||||
$ ABS(SX(I+2)) + ABS(SX(I+3)) +
|
||||
$ ABS(SX(I+4)) + ABS(SX(I+5))
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for increment not equal to 1
|
||||
*
|
||||
NINCX = N*INCX
|
||||
DO I = 1,NINCX,INCX
|
||||
STEMP = STEMP + ABS(SX(I))
|
||||
END DO
|
||||
END IF
|
||||
SASUM = STEMP
|
||||
RETURN
|
||||
*
|
||||
* End of SASUM
|
||||
*
|
||||
END
|
||||
*/
|
||||
float sasum ( int n, float x[], int incx ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
int i, j;
|
||||
float value;
|
||||
|
||||
value = 0.0;
|
||||
j = 0;
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
value = value + fabs ( x[j] );
|
||||
j = j + incx;
|
||||
}
|
||||
|
||||
return value;
|
||||
}
|
||||
130
src/saxpy.c
130
src/saxpy.c
|
|
@ -1,3 +1,6 @@
|
|||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b SAXPY
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -85,68 +88,65 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
REAL SA
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL SX(*),SY(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
INTEGER I,IX,IY,M,MP1
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MOD
|
||||
* ..
|
||||
IF (N.LE.0) RETURN
|
||||
IF (SA.EQ.0.0) RETURN
|
||||
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
|
||||
*
|
||||
* code for both increments equal to 1
|
||||
*
|
||||
*
|
||||
* clean-up loop
|
||||
*
|
||||
M = MOD(N,4)
|
||||
IF (M.NE.0) THEN
|
||||
DO I = 1,M
|
||||
SY(I) = SY(I) + SA*SX(I)
|
||||
END DO
|
||||
END IF
|
||||
IF (N.LT.4) RETURN
|
||||
MP1 = M + 1
|
||||
DO I = MP1,N,4
|
||||
SY(I) = SY(I) + SA*SX(I)
|
||||
SY(I+1) = SY(I+1) + SA*SX(I+1)
|
||||
SY(I+2) = SY(I+2) + SA*SX(I+2)
|
||||
SY(I+3) = SY(I+3) + SA*SX(I+3)
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for unequal increments or equal increments
|
||||
* not equal to 1
|
||||
*
|
||||
IX = 1
|
||||
IY = 1
|
||||
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
|
||||
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
|
||||
DO I = 1,N
|
||||
SY(IY) = SY(IY) + SA*SX(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
END DO
|
||||
END IF
|
||||
RETURN
|
||||
*
|
||||
* End of SAXPY
|
||||
*
|
||||
END
|
||||
*/
|
||||
void saxpy ( int n, float da, float dx[], int incx, float dy[], int incy ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
int i, ix, iy, m;
|
||||
|
||||
if ( n <= 0 ) {
|
||||
return;
|
||||
}
|
||||
|
||||
if ( da == 0.0 ) {
|
||||
return;
|
||||
}
|
||||
/*
|
||||
Code for unequal increments or equal increments
|
||||
not equal to 1.
|
||||
*/
|
||||
if ( incx != 1 || incy != 1 ) {
|
||||
if ( 0 <= incx ) {
|
||||
ix = 0;
|
||||
}
|
||||
else {
|
||||
ix = ( - n + 1 ) * incx;
|
||||
}
|
||||
|
||||
if ( 0 <= incy ) {
|
||||
iy = 0;
|
||||
}
|
||||
else {
|
||||
iy = ( - n + 1 ) * incy;
|
||||
}
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
dy[iy] = dy[iy] + da * dx[ix];
|
||||
ix = ix + incx;
|
||||
iy = iy + incy;
|
||||
}
|
||||
} else {
|
||||
// Code for both increments equal to 1.
|
||||
m = n % 4;
|
||||
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
dy[i] = dy[i] + da * dx[i];
|
||||
}
|
||||
|
||||
for ( i = m; i < n; i = i + 4 ) {
|
||||
dy[i ] = dy[i ] + da * dx[i ];
|
||||
dy[i+1] = dy[i+1] + da * dx[i+1];
|
||||
dy[i+2] = dy[i+2] + da * dx[i+2];
|
||||
dy[i+3] = dy[i+3] + da * dx[i+3];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
35
src/scabs1.c
35
src/scabs1.c
|
|
@ -1,3 +1,8 @@
|
|||
# include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b SCABS1
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -42,24 +47,12 @@
|
|||
*> \ingroup single_blas_level1
|
||||
*
|
||||
* =====================================================================
|
||||
REAL FUNCTION SCABS1(Z)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
COMPLEX Z
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS,AIMAG,REAL
|
||||
* ..
|
||||
SCABS1 = ABS(REAL(Z)) + ABS(AIMAG(Z))
|
||||
RETURN
|
||||
*
|
||||
* End of SCABS1
|
||||
*
|
||||
END
|
||||
*/
|
||||
float cabs1 ( float complex z ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
return fabs ( creal ( z ) ) + fabs ( cimag ( z ) );
|
||||
}
|
||||
|
|
|
|||
90
src/scasum.c
90
src/scasum.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b SCASUM
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -68,50 +73,41 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
REAL FUNCTION SCASUM(N,CX,INCX)
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX CX(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
REAL STEMP
|
||||
INTEGER I,NINCX
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS,AIMAG,REAL
|
||||
* ..
|
||||
SCASUM = 0.0e0
|
||||
STEMP = 0.0e0
|
||||
IF (N.LE.0 .OR. INCX.LE.0) RETURN
|
||||
IF (INCX.EQ.1) THEN
|
||||
*
|
||||
* code for increment equal to 1
|
||||
*
|
||||
DO I = 1,N
|
||||
STEMP = STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I)))
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* code for increment not equal to 1
|
||||
*
|
||||
NINCX = N*INCX
|
||||
DO I = 1,NINCX,INCX
|
||||
STEMP = STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I)))
|
||||
END DO
|
||||
END IF
|
||||
SCASUM = STEMP
|
||||
RETURN
|
||||
*
|
||||
* End of SCASUM
|
||||
*
|
||||
END
|
||||
*/
|
||||
float scasum ( int n, float complex x[], int incx ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
int i, ix;
|
||||
float value;
|
||||
|
||||
value = 0.0;
|
||||
|
||||
if ( n <= 0 || incx <= 0 ) {
|
||||
return value;
|
||||
}
|
||||
|
||||
if ( incx == 1 ) {
|
||||
for ( i = 0; i < n; i++ )
|
||||
{
|
||||
value = value + fabs ( creal ( x[i] ) )
|
||||
+ fabs ( cimag ( x[i] ) );
|
||||
}
|
||||
}
|
||||
else {
|
||||
ix = 0;
|
||||
for ( i = 0; i < n; i++ )
|
||||
{
|
||||
value = value + fabs ( creal ( x[ix] ) )
|
||||
+ fabs ( cimag ( x[ix] ) );
|
||||
ix = ix + incx;
|
||||
}
|
||||
}
|
||||
|
||||
return value;
|
||||
}
|
||||
351
src/scnrm2.c
351
src/scnrm2.c
|
|
@ -1,209 +1,144 @@
|
|||
!> \brief \b SCNRM2
|
||||
!
|
||||
! =========== DOCUMENTATION ===========
|
||||
!
|
||||
! Online html documentation available at
|
||||
! http://www.netlib.org/lapack/explore-html/
|
||||
!
|
||||
! Definition:
|
||||
! ===========
|
||||
!
|
||||
! REAL FUNCTION SCNRM2(N,X,INCX)
|
||||
!
|
||||
! .. Scalar Arguments ..
|
||||
! INTEGER INCX,N
|
||||
! ..
|
||||
! .. Array Arguments ..
|
||||
! COMPLEX X(*)
|
||||
! ..
|
||||
!
|
||||
!
|
||||
!> \par Purpose:
|
||||
! =============
|
||||
!>
|
||||
!> \verbatim
|
||||
!>
|
||||
!> SCNRM2 returns the euclidean norm of a vector via the function
|
||||
!> name, so that
|
||||
!>
|
||||
!> SCNRM2 := sqrt( x**H*x )
|
||||
!> \endverbatim
|
||||
!
|
||||
! Arguments:
|
||||
! ==========
|
||||
!
|
||||
!> \param[in] N
|
||||
!> \verbatim
|
||||
!> N is INTEGER
|
||||
!> number of elements in input vector(s)
|
||||
!> \endverbatim
|
||||
!>
|
||||
!> \param[in] X
|
||||
!> \verbatim
|
||||
!> X is COMPLEX array, dimension (N)
|
||||
!> complex vector with N elements
|
||||
!> \endverbatim
|
||||
!>
|
||||
!> \param[in] INCX
|
||||
!> \verbatim
|
||||
!> INCX is INTEGER, storage spacing between elements of X
|
||||
!> If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
|
||||
!> If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
|
||||
!> If INCX = 0, x isn't a vector so there is no need to call
|
||||
!> this subroutine. If you call it anyway, it will count x(1)
|
||||
!> in the vector norm N times.
|
||||
!> \endverbatim
|
||||
!
|
||||
! Authors:
|
||||
! ========
|
||||
!
|
||||
!> \author Edward Anderson, Lockheed Martin
|
||||
!
|
||||
!> \date August 2016
|
||||
!
|
||||
!> \ingroup single_blas_level1
|
||||
!
|
||||
!> \par Contributors:
|
||||
! ==================
|
||||
!>
|
||||
!> Weslley Pereira, University of Colorado Denver, USA
|
||||
!
|
||||
!> \par Further Details:
|
||||
! =====================
|
||||
!>
|
||||
!> \verbatim
|
||||
!>
|
||||
!> Anderson E. (2017)
|
||||
!> Algorithm 978: Safe Scaling in the Level 1 BLAS
|
||||
!> ACM Trans Math Softw 44:1--28
|
||||
!> https://doi.org/10.1145/3061665
|
||||
!>
|
||||
!> Blue, James L. (1978)
|
||||
!> A Portable Fortran Program to Find the Euclidean Norm of a Vector
|
||||
!> ACM Trans Math Softw 4:15--23
|
||||
!> https://doi.org/10.1145/355769.355771
|
||||
!>
|
||||
!> \endverbatim
|
||||
!>
|
||||
! =====================================================================
|
||||
function SCNRM2( n, x, incx )
|
||||
integer, parameter :: wp = kind(1.e0)
|
||||
real(wp) :: SCNRM2
|
||||
!
|
||||
! -- Reference BLAS level1 routine (version 3.9.1) --
|
||||
! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
! March 2021
|
||||
!
|
||||
! .. Constants ..
|
||||
real(wp), parameter :: zero = 0.0_wp
|
||||
real(wp), parameter :: one = 1.0_wp
|
||||
real(wp), parameter :: maxN = huge(0.0_wp)
|
||||
! ..
|
||||
! .. Blue's scaling constants ..
|
||||
real(wp), parameter :: tsml = real(radix(0._wp), wp)**ceiling( &
|
||||
(minexponent(0._wp) - 1) * 0.5_wp)
|
||||
real(wp), parameter :: tbig = real(radix(0._wp), wp)**floor( &
|
||||
(maxexponent(0._wp) - digits(0._wp) + 1) * 0.5_wp)
|
||||
real(wp), parameter :: ssml = real(radix(0._wp), wp)**( - floor( &
|
||||
(minexponent(0._wp) - digits(0._wp)) * 0.5_wp))
|
||||
real(wp), parameter :: sbig = real(radix(0._wp), wp)**( - ceiling( &
|
||||
(maxexponent(0._wp) + digits(0._wp) - 1) * 0.5_wp))
|
||||
! ..
|
||||
! .. Scalar Arguments ..
|
||||
integer :: incx, n
|
||||
! ..
|
||||
! .. Array Arguments ..
|
||||
complex(wp) :: x(*)
|
||||
! ..
|
||||
! .. Local Scalars ..
|
||||
integer :: i, ix
|
||||
logical :: notbig
|
||||
real(wp) :: abig, amed, asml, ax, scl, sumsq, ymax, ymin
|
||||
!
|
||||
! Quick return if possible
|
||||
!
|
||||
SCNRM2 = zero
|
||||
if( n <= 0 ) return
|
||||
!
|
||||
scl = one
|
||||
sumsq = zero
|
||||
!
|
||||
! Compute the sum of squares in 3 accumulators:
|
||||
! abig -- sums of squares scaled down to avoid overflow
|
||||
! asml -- sums of squares scaled up to avoid underflow
|
||||
! amed -- sums of squares that do not require scaling
|
||||
! The thresholds and multipliers are
|
||||
! tbig -- values bigger than this are scaled down by sbig
|
||||
! tsml -- values smaller than this are scaled up by ssml
|
||||
!
|
||||
notbig = .true.
|
||||
asml = zero
|
||||
amed = zero
|
||||
abig = zero
|
||||
ix = 1
|
||||
if( incx < 0 ) ix = 1 - (n-1)*incx
|
||||
do i = 1, n
|
||||
ax = abs(real(x(ix)))
|
||||
if (ax > tbig) then
|
||||
abig = abig + (ax*sbig)**2
|
||||
notbig = .false.
|
||||
else if (ax < tsml) then
|
||||
if (notbig) asml = asml + (ax*ssml)**2
|
||||
#include <complex.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b SCNRM2
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
*
|
||||
* Online html documentation available at
|
||||
* http://www.netlib.org/lapack/explore-html/
|
||||
*
|
||||
* Definition:
|
||||
* ===========
|
||||
*
|
||||
* REAL FUNCTION SCNRM2(N,X,INCX)
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
* INTEGER INCX,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
* COMPLEX X(*)
|
||||
* ..
|
||||
*
|
||||
*
|
||||
*> \par Purpose:
|
||||
* =============
|
||||
*>
|
||||
*> \verbatim
|
||||
*>
|
||||
*> SCNRM2 returns the euclidean norm of a vector via the function
|
||||
*> name, so that
|
||||
*>
|
||||
*> SCNRM2 := sqrt( x**H*x )
|
||||
*> \endverbatim
|
||||
*
|
||||
* Arguments:
|
||||
* ==========
|
||||
*
|
||||
*> \param[in] N
|
||||
*> \verbatim
|
||||
*> N is INTEGER
|
||||
*> number of elements in input vector(s)
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] X
|
||||
*> \verbatim
|
||||
*> X is COMPLEX array, dimension (N)
|
||||
*> complex vector with N elements
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] INCX
|
||||
*> \verbatim
|
||||
*> INCX is INTEGER, storage spacing between elements of X
|
||||
*> If INCX > 0, X(1+(i-1)*INCX) = x(i) for 1 <= i <= n
|
||||
*> If INCX < 0, X(1-(n-i)*INCX) = x(i) for 1 <= i <= n
|
||||
*> If INCX = 0, x isn't a vector so there is no need to call
|
||||
*> this subroutine. If you call it anyway, it will count x(1)
|
||||
*> in the vector norm N times.
|
||||
*> \endverbatim
|
||||
*
|
||||
* Authors:
|
||||
* ========
|
||||
*
|
||||
*> \author Edward Anderson, Lockheed Martin
|
||||
*
|
||||
*> \date August 2016
|
||||
*
|
||||
*> \ingroup single_blas_level1
|
||||
*
|
||||
*> \par Contributors:
|
||||
* ==================
|
||||
*>
|
||||
*> Weslley Pereira, University of Colorado Denver, USA
|
||||
*
|
||||
*> \par Further Details:
|
||||
* =====================
|
||||
*>
|
||||
*> \verbatim
|
||||
*>
|
||||
*> Anderson E. (2017)
|
||||
*> Algorithm 978: Safe Scaling in the Level 1 BLAS
|
||||
*> ACM Trans Math Softw 44:1--28
|
||||
*> https://doi.org/10.1145/3061665
|
||||
*>
|
||||
*> Blue, James L. (1978)
|
||||
*> A Portable Fortran Program to Find the Euclidean Norm of a Vector
|
||||
*> ACM Trans Math Softw 4:15--23
|
||||
*> https://doi.org/10.1145/355769.355771
|
||||
*>
|
||||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
*/
|
||||
float scnrm2 ( int n, float complex x[], int incx ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
int i, ix;
|
||||
float scale, ssq, temp, value;
|
||||
|
||||
if ( n < 1 || incx < 1 ) {
|
||||
value = 0.0;
|
||||
} else {
|
||||
scale = 0.0;
|
||||
ssq = 1.0;
|
||||
ix = 0;
|
||||
|
||||
for ( i = 0; i < n; i++ ) {
|
||||
if ( creal ( x[ix] ) != 0.0 ) {
|
||||
temp = fabs ( creal ( x[ix] ) );
|
||||
if ( scale < temp ) {
|
||||
ssq = 1.0 + ssq * pow ( scale / temp, 2 );
|
||||
scale = temp;
|
||||
} else {
|
||||
ssq = ssq + pow ( temp / scale, 2 );
|
||||
}
|
||||
}
|
||||
|
||||
if ( cimag ( x[ix] ) != 0.0 ) {
|
||||
temp = fabs ( cimag ( x[ix] ) );
|
||||
if ( scale < temp )
|
||||
{
|
||||
ssq = 1.0 + ssq * pow ( scale / temp, 2 );
|
||||
scale = temp;
|
||||
}
|
||||
else
|
||||
amed = amed + ax**2
|
||||
end if
|
||||
ax = abs(aimag(x(ix)))
|
||||
if (ax > tbig) then
|
||||
abig = abig + (ax*sbig)**2
|
||||
notbig = .false.
|
||||
else if (ax < tsml) then
|
||||
if (notbig) asml = asml + (ax*ssml)**2
|
||||
else
|
||||
amed = amed + ax**2
|
||||
end if
|
||||
ix = ix + incx
|
||||
end do
|
||||
!
|
||||
! Combine abig and amed or amed and asml if more than one
|
||||
! accumulator was used.
|
||||
!
|
||||
if (abig > zero) then
|
||||
!
|
||||
! Combine abig and amed if abig > 0.
|
||||
!
|
||||
if ( (amed > zero) .or. (amed > maxN) .or. (amed /= amed) ) then
|
||||
abig = abig + (amed*sbig)*sbig
|
||||
end if
|
||||
scl = one / sbig
|
||||
sumsq = abig
|
||||
else if (asml > zero) then
|
||||
!
|
||||
! Combine amed and asml if asml > 0.
|
||||
!
|
||||
if ( (amed > zero) .or. (amed > maxN) .or. (amed /= amed) ) then
|
||||
amed = sqrt(amed)
|
||||
asml = sqrt(asml) / ssml
|
||||
if (asml > amed) then
|
||||
ymin = amed
|
||||
ymax = asml
|
||||
else
|
||||
ymin = asml
|
||||
ymax = amed
|
||||
end if
|
||||
scl = one
|
||||
sumsq = ymax**2*( one + (ymin/ymax)**2 )
|
||||
else
|
||||
scl = one / ssml
|
||||
sumsq = asml
|
||||
end if
|
||||
else
|
||||
!
|
||||
! Otherwise all values are mid-range
|
||||
!
|
||||
scl = one
|
||||
sumsq = amed
|
||||
end if
|
||||
SCNRM2 = scl*sqrt( sumsq )
|
||||
return
|
||||
end function
|
||||
{
|
||||
ssq = ssq + pow ( temp / scale, 2 );
|
||||
}
|
||||
}
|
||||
ix = ix + incx;
|
||||
}
|
||||
value = scale * sqrt ( ssq );
|
||||
}
|
||||
|
||||
return value;
|
||||
|
||||
}
|
||||
337
src/sgemv.c
337
src/sgemv.c
|
|
@ -1,3 +1,6 @@
|
|||
# include "blas0.h"
|
||||
|
||||
/*
|
||||
*> \brief \b SGEMV
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -152,176 +155,164 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE SGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
|
||||
*
|
||||
* -- Reference BLAS level2 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
REAL ALPHA,BETA
|
||||
INTEGER INCX,INCY,LDA,M,N
|
||||
CHARACTER TRANS
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL A(LDA,*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
REAL ONE,ZERO
|
||||
PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
REAL TEMP
|
||||
INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MAX
|
||||
* ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
|
||||
+ .NOT.LSAME(TRANS,'C')) THEN
|
||||
INFO = 1
|
||||
ELSE IF (M.LT.0) THEN
|
||||
INFO = 2
|
||||
ELSE IF (N.LT.0) THEN
|
||||
INFO = 3
|
||||
ELSE IF (LDA.LT.MAX(1,M)) THEN
|
||||
INFO = 6
|
||||
ELSE IF (INCX.EQ.0) THEN
|
||||
INFO = 8
|
||||
ELSE IF (INCY.EQ.0) THEN
|
||||
INFO = 11
|
||||
END IF
|
||||
IF (INFO.NE.0) THEN
|
||||
CALL XERBLA('SGEMV ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
|
||||
+ ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
|
||||
*
|
||||
* Set LENX and LENY, the lengths of the vectors x and y, and set
|
||||
* up the start points in X and Y.
|
||||
*
|
||||
IF (LSAME(TRANS,'N')) THEN
|
||||
LENX = N
|
||||
LENY = M
|
||||
ELSE
|
||||
LENX = M
|
||||
LENY = N
|
||||
END IF
|
||||
IF (INCX.GT.0) THEN
|
||||
KX = 1
|
||||
ELSE
|
||||
KX = 1 - (LENX-1)*INCX
|
||||
END IF
|
||||
IF (INCY.GT.0) THEN
|
||||
KY = 1
|
||||
ELSE
|
||||
KY = 1 - (LENY-1)*INCY
|
||||
END IF
|
||||
*
|
||||
* Start the operations. In this version the elements of A are
|
||||
* accessed sequentially with one pass through A.
|
||||
*
|
||||
* First form y := beta*y.
|
||||
*
|
||||
IF (BETA.NE.ONE) THEN
|
||||
IF (INCY.EQ.1) THEN
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
DO 10 I = 1,LENY
|
||||
Y(I) = ZERO
|
||||
10 CONTINUE
|
||||
ELSE
|
||||
DO 20 I = 1,LENY
|
||||
Y(I) = BETA*Y(I)
|
||||
20 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IY = KY
|
||||
IF (BETA.EQ.ZERO) THEN
|
||||
DO 30 I = 1,LENY
|
||||
Y(IY) = ZERO
|
||||
IY = IY + INCY
|
||||
30 CONTINUE
|
||||
ELSE
|
||||
DO 40 I = 1,LENY
|
||||
Y(IY) = BETA*Y(IY)
|
||||
IY = IY + INCY
|
||||
40 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
IF (ALPHA.EQ.ZERO) RETURN
|
||||
IF (LSAME(TRANS,'N')) THEN
|
||||
*
|
||||
* Form y := alpha*A*x + y.
|
||||
*
|
||||
JX = KX
|
||||
IF (INCY.EQ.1) THEN
|
||||
DO 60 J = 1,N
|
||||
TEMP = ALPHA*X(JX)
|
||||
DO 50 I = 1,M
|
||||
Y(I) = Y(I) + TEMP*A(I,J)
|
||||
50 CONTINUE
|
||||
JX = JX + INCX
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
DO 80 J = 1,N
|
||||
TEMP = ALPHA*X(JX)
|
||||
IY = KY
|
||||
DO 70 I = 1,M
|
||||
Y(IY) = Y(IY) + TEMP*A(I,J)
|
||||
IY = IY + INCY
|
||||
70 CONTINUE
|
||||
JX = JX + INCX
|
||||
80 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
*
|
||||
* Form y := alpha*A**T*x + y.
|
||||
*
|
||||
JY = KY
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 100 J = 1,N
|
||||
TEMP = ZERO
|
||||
DO 90 I = 1,M
|
||||
TEMP = TEMP + A(I,J)*X(I)
|
||||
90 CONTINUE
|
||||
Y(JY) = Y(JY) + ALPHA*TEMP
|
||||
JY = JY + INCY
|
||||
100 CONTINUE
|
||||
ELSE
|
||||
DO 120 J = 1,N
|
||||
TEMP = ZERO
|
||||
IX = KX
|
||||
DO 110 I = 1,M
|
||||
TEMP = TEMP + A(I,J)*X(IX)
|
||||
IX = IX + INCX
|
||||
110 CONTINUE
|
||||
Y(JY) = Y(JY) + ALPHA*TEMP
|
||||
JY = JY + INCY
|
||||
120 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of SGEMV
|
||||
*
|
||||
END
|
||||
*/
|
||||
void sgemv ( char trans, int m, int n, float alpha, float a[], int lda,
|
||||
float x[], int incx, float beta, float y[], int incy ) {
|
||||
|
||||
// -- Reference BLAS level2 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
float temp;
|
||||
int i, info, ix, iy, j, jx, jy, kx, ky, lenx, leny;
|
||||
|
||||
// Test the input parameters.
|
||||
|
||||
info = 0;
|
||||
if ( ! lsame ( trans, 'N' ) &&
|
||||
! lsame ( trans, 'T' ) &&
|
||||
! lsame ( trans, 'C' ) ) {
|
||||
info = 1;
|
||||
} else if ( m < 0 ) {
|
||||
info = 2;
|
||||
} else if ( n < 0 ) {
|
||||
info = 3;
|
||||
} else if ( lda < i4_max ( 1, m ) ) {
|
||||
info = 6;
|
||||
} else if ( incx == 0 ) {
|
||||
info = 8;
|
||||
} else if ( incy == 0 ) {
|
||||
info = 11;
|
||||
}
|
||||
|
||||
if ( info != 0 ) {
|
||||
xerbla ( "SGEMV", info );
|
||||
return;
|
||||
}
|
||||
|
||||
// Quick return if possible.
|
||||
|
||||
if ( ( m == 0 ) ||
|
||||
( n == 0 ) ||
|
||||
( ( alpha == 0.0 ) && ( beta == 1.0 ) ) ) {
|
||||
return;
|
||||
}
|
||||
/*
|
||||
Set LENX and LENY, the lengths of the vectors x and y, and set
|
||||
up the start points in X and Y.
|
||||
*/
|
||||
if ( lsame ( trans, 'N' ) ) {
|
||||
lenx = n;
|
||||
leny = m;
|
||||
} else {
|
||||
lenx = m;
|
||||
leny = n;
|
||||
}
|
||||
|
||||
if ( 0 < incx ) {
|
||||
kx = 0;
|
||||
} else {
|
||||
kx = 0 - ( lenx - 1 ) * incx;
|
||||
}
|
||||
|
||||
if ( 0 < incy ) {
|
||||
ky = 0;
|
||||
} else {
|
||||
ky = 0 - ( leny - 1 ) * incy;
|
||||
}
|
||||
/*
|
||||
Start the operations. In this version the elements of A are
|
||||
accessed sequentially with one pass through A.
|
||||
|
||||
First form y := beta*y.
|
||||
*/
|
||||
if ( beta != 1.0 ) {
|
||||
if ( incy == 1 ) {
|
||||
if ( beta == 0.0 ) {
|
||||
for ( i = 0; i < leny; i++ ) {
|
||||
y[i] = 0.0;
|
||||
}
|
||||
} else {
|
||||
for ( i = 0; i < leny; i++ ) {
|
||||
y[i] = beta * y[i];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
iy = ky;
|
||||
if ( beta == 0.0 ) {
|
||||
for ( i = 0; i < leny; i++ ) {
|
||||
y[iy] = 0.0;
|
||||
iy = iy + incy;
|
||||
}
|
||||
} else {
|
||||
for ( i = 0; i < leny; i++ ) {
|
||||
y[iy] = beta * y[iy];
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if ( alpha == 0.0 ) {
|
||||
return;
|
||||
}
|
||||
|
||||
// Form y := alpha*A*x + y.
|
||||
|
||||
if ( lsame ( trans, 'N' ) ) {
|
||||
jx = kx;
|
||||
if ( incy == 1 ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( x[jx] != 0.0 ) {
|
||||
temp = alpha * x[jx];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
y[i] = y[i] + temp * a[i+j*lda];
|
||||
}
|
||||
}
|
||||
jx = jx + incx;
|
||||
}
|
||||
} else {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( x[jx] != 0.0 ) {
|
||||
temp = alpha * x[jx];
|
||||
iy = ky;
|
||||
for ( i = 0; i < m; i++ )
|
||||
{
|
||||
y[iy] = y[iy] + temp * a[i+j*lda];
|
||||
iy = iy + incy;
|
||||
}
|
||||
}
|
||||
jx = jx + incx;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// Form y := alpha*A'*x + y.
|
||||
jy = ky;
|
||||
if ( incx == 1 ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
temp = 0.0;
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = temp + a[i+j*lda] * x[i];
|
||||
}
|
||||
y[jy] = y[jy] + alpha * temp;
|
||||
jy = jy + incy;
|
||||
}
|
||||
} else {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
temp = 0.0;
|
||||
ix = kx;
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = temp + a[i+j*lda] * x[ix];
|
||||
ix = ix + incx;
|
||||
}
|
||||
y[jy] = y[jy] + alpha * temp;
|
||||
jy = jy + incy;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
181
src/sger.c
181
src/sger.c
|
|
@ -1,3 +1,6 @@
|
|||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b SGER
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -126,99 +129,85 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE SGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
|
||||
*
|
||||
* -- Reference BLAS level2 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
REAL ALPHA
|
||||
INTEGER INCX,INCY,LDA,M,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL A(LDA,*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
REAL ZERO
|
||||
PARAMETER (ZERO=0.0E+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
REAL TEMP
|
||||
INTEGER I,INFO,IX,J,JY,KX
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MAX
|
||||
* ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
IF (M.LT.0) THEN
|
||||
INFO = 1
|
||||
ELSE IF (N.LT.0) THEN
|
||||
INFO = 2
|
||||
ELSE IF (INCX.EQ.0) THEN
|
||||
INFO = 5
|
||||
ELSE IF (INCY.EQ.0) THEN
|
||||
INFO = 7
|
||||
ELSE IF (LDA.LT.MAX(1,M)) THEN
|
||||
INFO = 9
|
||||
END IF
|
||||
IF (INFO.NE.0) THEN
|
||||
CALL XERBLA('SGER ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF ((M.EQ.0) .OR. (N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
|
||||
*
|
||||
* Start the operations. In this version the elements of A are
|
||||
* accessed sequentially with one pass through A.
|
||||
*
|
||||
IF (INCY.GT.0) THEN
|
||||
JY = 1
|
||||
ELSE
|
||||
JY = 1 - (N-1)*INCY
|
||||
END IF
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 20 J = 1,N
|
||||
IF (Y(JY).NE.ZERO) THEN
|
||||
TEMP = ALPHA*Y(JY)
|
||||
DO 10 I = 1,M
|
||||
A(I,J) = A(I,J) + X(I)*TEMP
|
||||
10 CONTINUE
|
||||
END IF
|
||||
JY = JY + INCY
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
IF (INCX.GT.0) THEN
|
||||
KX = 1
|
||||
ELSE
|
||||
KX = 1 - (M-1)*INCX
|
||||
END IF
|
||||
DO 40 J = 1,N
|
||||
IF (Y(JY).NE.ZERO) THEN
|
||||
TEMP = ALPHA*Y(JY)
|
||||
IX = KX
|
||||
DO 30 I = 1,M
|
||||
A(I,J) = A(I,J) + X(IX)*TEMP
|
||||
IX = IX + INCX
|
||||
30 CONTINUE
|
||||
END IF
|
||||
JY = JY + INCY
|
||||
40 CONTINUE
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of SGER
|
||||
*
|
||||
END
|
||||
*/
|
||||
void sger ( int m, int n, float alpha, float x[], int incx, float y[],
|
||||
int incy, float a[], int lda ) {
|
||||
|
||||
// -- Reference BLAS level2 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
float temp;
|
||||
int i, info, ix, j, jy, kx;
|
||||
|
||||
// Test the input parameters.
|
||||
|
||||
info = 0;
|
||||
if ( m < 0 ) {
|
||||
info = 1;
|
||||
} else if ( n < 0 ) {
|
||||
info = 2;
|
||||
} else if ( incx == 0 ) {
|
||||
info = 5;
|
||||
} else if ( incy == 0 ) {
|
||||
info = 7;
|
||||
} else if ( lda < i4_max ( 1, m ) ) {
|
||||
info = 9;
|
||||
}
|
||||
|
||||
if ( info != 0 ) {
|
||||
xerbla ( "SGER", info );
|
||||
return;
|
||||
}
|
||||
|
||||
// Quick return if possible.
|
||||
|
||||
if ( m == 0 || n == 0 || alpha == 0.0 ) {
|
||||
return;
|
||||
}
|
||||
/*
|
||||
Start the operations. In this version the elements of A are
|
||||
accessed sequentially with one pass through A.
|
||||
*/
|
||||
if ( 0 < incy ) {
|
||||
jy = 0;
|
||||
} else {
|
||||
jy = 0 - ( n - 1 ) * incy;
|
||||
}
|
||||
|
||||
if ( incx == 1 ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( y[jy] != 0.0 ) {
|
||||
temp = alpha * y[jy];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
a[i+j*lda] = a[i+j*lda] + x[i] * temp;
|
||||
}
|
||||
}
|
||||
jy = jy + incy;
|
||||
}
|
||||
} else {
|
||||
if ( 0 < incx ) {
|
||||
kx = 0;
|
||||
}
|
||||
else {
|
||||
kx = 0 - ( m - 1 ) * incx;
|
||||
}
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( y[jy] != 0.0 ) {
|
||||
temp = alpha * y[jy];
|
||||
ix = kx;
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
a[i+j*lda] = a[i+j*lda] + x[ix] * temp;
|
||||
ix = ix + incx;
|
||||
}
|
||||
}
|
||||
jy = jy + incy;
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
}
|
||||
384
src/strmv.c
384
src/strmv.c
|
|
@ -1,3 +1,6 @@
|
|||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b STRMV
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -143,197 +146,190 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE STRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
|
||||
*
|
||||
* -- Reference BLAS level2 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,LDA,N
|
||||
CHARACTER DIAG,TRANS,UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL A(LDA,*),X(*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
REAL ZERO
|
||||
PARAMETER (ZERO=0.0E+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
REAL TEMP
|
||||
INTEGER I,INFO,IX,J,JX,KX
|
||||
LOGICAL NOUNIT
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MAX
|
||||
* ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
|
||||
INFO = 1
|
||||
ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
|
||||
+ .NOT.LSAME(TRANS,'C')) THEN
|
||||
INFO = 2
|
||||
ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
|
||||
INFO = 3
|
||||
ELSE IF (N.LT.0) THEN
|
||||
INFO = 4
|
||||
ELSE IF (LDA.LT.MAX(1,N)) THEN
|
||||
INFO = 6
|
||||
ELSE IF (INCX.EQ.0) THEN
|
||||
INFO = 8
|
||||
END IF
|
||||
IF (INFO.NE.0) THEN
|
||||
CALL XERBLA('STRMV ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF (N.EQ.0) RETURN
|
||||
*
|
||||
NOUNIT = LSAME(DIAG,'N')
|
||||
*
|
||||
* Set up the start point in X if the increment is not unity. This
|
||||
* will be ( N - 1 )*INCX too small for descending loops.
|
||||
*
|
||||
IF (INCX.LE.0) THEN
|
||||
KX = 1 - (N-1)*INCX
|
||||
ELSE IF (INCX.NE.1) THEN
|
||||
KX = 1
|
||||
END IF
|
||||
*
|
||||
* Start the operations. In this version the elements of A are
|
||||
* accessed sequentially with one pass through A.
|
||||
*
|
||||
IF (LSAME(TRANS,'N')) THEN
|
||||
*
|
||||
* Form x := A*x.
|
||||
*
|
||||
IF (LSAME(UPLO,'U')) THEN
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 20 J = 1,N
|
||||
IF (X(J).NE.ZERO) THEN
|
||||
TEMP = X(J)
|
||||
DO 10 I = 1,J - 1
|
||||
X(I) = X(I) + TEMP*A(I,J)
|
||||
10 CONTINUE
|
||||
IF (NOUNIT) X(J) = X(J)*A(J,J)
|
||||
END IF
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 40 J = 1,N
|
||||
IF (X(JX).NE.ZERO) THEN
|
||||
TEMP = X(JX)
|
||||
IX = KX
|
||||
DO 30 I = 1,J - 1
|
||||
X(IX) = X(IX) + TEMP*A(I,J)
|
||||
IX = IX + INCX
|
||||
30 CONTINUE
|
||||
IF (NOUNIT) X(JX) = X(JX)*A(J,J)
|
||||
END IF
|
||||
JX = JX + INCX
|
||||
40 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 60 J = N,1,-1
|
||||
IF (X(J).NE.ZERO) THEN
|
||||
TEMP = X(J)
|
||||
DO 50 I = N,J + 1,-1
|
||||
X(I) = X(I) + TEMP*A(I,J)
|
||||
50 CONTINUE
|
||||
IF (NOUNIT) X(J) = X(J)*A(J,J)
|
||||
END IF
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
KX = KX + (N-1)*INCX
|
||||
JX = KX
|
||||
DO 80 J = N,1,-1
|
||||
IF (X(JX).NE.ZERO) THEN
|
||||
TEMP = X(JX)
|
||||
IX = KX
|
||||
DO 70 I = N,J + 1,-1
|
||||
X(IX) = X(IX) + TEMP*A(I,J)
|
||||
IX = IX - INCX
|
||||
70 CONTINUE
|
||||
IF (NOUNIT) X(JX) = X(JX)*A(J,J)
|
||||
END IF
|
||||
JX = JX - INCX
|
||||
80 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
ELSE
|
||||
*
|
||||
* Form x := A**T*x.
|
||||
*
|
||||
IF (LSAME(UPLO,'U')) THEN
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 100 J = N,1,-1
|
||||
TEMP = X(J)
|
||||
IF (NOUNIT) TEMP = TEMP*A(J,J)
|
||||
DO 90 I = J - 1,1,-1
|
||||
TEMP = TEMP + A(I,J)*X(I)
|
||||
90 CONTINUE
|
||||
X(J) = TEMP
|
||||
100 CONTINUE
|
||||
ELSE
|
||||
JX = KX + (N-1)*INCX
|
||||
DO 120 J = N,1,-1
|
||||
TEMP = X(JX)
|
||||
IX = JX
|
||||
IF (NOUNIT) TEMP = TEMP*A(J,J)
|
||||
DO 110 I = J - 1,1,-1
|
||||
IX = IX - INCX
|
||||
TEMP = TEMP + A(I,J)*X(IX)
|
||||
110 CONTINUE
|
||||
X(JX) = TEMP
|
||||
JX = JX - INCX
|
||||
120 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 140 J = 1,N
|
||||
TEMP = X(J)
|
||||
IF (NOUNIT) TEMP = TEMP*A(J,J)
|
||||
DO 130 I = J + 1,N
|
||||
TEMP = TEMP + A(I,J)*X(I)
|
||||
130 CONTINUE
|
||||
X(J) = TEMP
|
||||
140 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 160 J = 1,N
|
||||
TEMP = X(JX)
|
||||
IX = JX
|
||||
IF (NOUNIT) TEMP = TEMP*A(J,J)
|
||||
DO 150 I = J + 1,N
|
||||
IX = IX + INCX
|
||||
TEMP = TEMP + A(I,J)*X(IX)
|
||||
150 CONTINUE
|
||||
X(JX) = TEMP
|
||||
JX = JX + INCX
|
||||
160 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of STRMV
|
||||
*
|
||||
END
|
||||
*/
|
||||
void strmv(char uplo, char trans, char diag, int n, float a[], int lda,
|
||||
float x[], int incx) {
|
||||
// -- Reference BLAS level2 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
int i, info, ix, j, jx, kx, nounit;
|
||||
float temp;
|
||||
|
||||
// Test the input parameters.
|
||||
|
||||
info = 0;
|
||||
if ( ! lsame ( uplo, 'U' ) && ! lsame ( uplo, 'L' ) ) {
|
||||
info = 1;
|
||||
} else if ( ! lsame ( trans, 'N' ) && ! lsame ( trans, 'T' ) &&
|
||||
! lsame ( trans, 'C' ) ) {
|
||||
info = 2;
|
||||
} else if ( ! lsame ( diag, 'U' ) && ! lsame ( diag, 'N' ) ) {
|
||||
info = 3;
|
||||
} else if ( n < 0 ) {
|
||||
info = 4;
|
||||
} else if ( lda < i4_max ( 1, n ) ) {
|
||||
info = 6;
|
||||
} else if ( incx == 0 ) {
|
||||
info = 8;
|
||||
}
|
||||
|
||||
if ( info != 0 ) {
|
||||
xerbla ( "STRMV", info );
|
||||
return;
|
||||
}
|
||||
|
||||
// Quick return if possible.
|
||||
|
||||
if ( n == 0 ) {
|
||||
return;
|
||||
}
|
||||
|
||||
nounit = lsame ( diag, 'N' );
|
||||
/*
|
||||
Set up the start point in X if the increment is not unity. This
|
||||
will be ( N - 1 ) * INCX too small for descending loops.
|
||||
*/
|
||||
if ( incx <= 0 ) {
|
||||
kx = 0 - ( n - 1 ) * incx;
|
||||
} else if ( incx != 1 ) {
|
||||
kx = 0;
|
||||
}
|
||||
/*
|
||||
Start the operations. In this version the elements of A are
|
||||
accessed sequentially with one pass through A.
|
||||
*/
|
||||
if ( lsame ( trans, 'N' ) ) {
|
||||
|
||||
// Form x := A*x.
|
||||
|
||||
if ( lsame ( uplo, 'U' ) ) {
|
||||
if ( incx == 1 ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( x[j] != 0.0 ) {
|
||||
temp = x[j];
|
||||
for ( i = 0; i < j; i++ ) {
|
||||
x[i] = x[i] + temp * a[i+j*lda];
|
||||
}
|
||||
if ( nounit ) {
|
||||
x[j] = x[j] * a[j+j*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
jx = kx;
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( x[jx] != 0.0 ) {
|
||||
temp = x[jx];
|
||||
ix = kx;
|
||||
for ( i = 0; i < j; i++ ) {
|
||||
x[ix] = x[ix] + temp * a[i+j*lda];
|
||||
ix = ix + incx;
|
||||
}
|
||||
if ( nounit ) {
|
||||
x[jx] = x[jx] * a[j+j*lda];
|
||||
}
|
||||
}
|
||||
jx = jx + incx;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if ( incx == 1 ) {
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
if ( x[j] != 0.0 ) {
|
||||
temp = x[j];
|
||||
for ( i = n - 1; j < i; i-- ) {
|
||||
x[i] = x[i] + temp * a[i+j*lda];
|
||||
}
|
||||
if ( nounit ) {
|
||||
x[j] = x[j] * a[j+j*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
kx = kx + ( n - 1 ) * incx;
|
||||
jx = kx;
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
if ( x[jx] != 0.0 ) {
|
||||
temp = x[jx];
|
||||
ix = kx;
|
||||
for ( i = n - 1; j < i; i-- ) {
|
||||
x[ix] = x[ix] + temp * a[i+j*lda];
|
||||
ix = ix - incx;
|
||||
}
|
||||
if ( nounit ) {
|
||||
x[jx] = x[jx] * a[j+j*lda];
|
||||
}
|
||||
}
|
||||
jx = jx - incx;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// Form x := A'*x.
|
||||
if ( lsame ( uplo, 'U' ) ) {
|
||||
if ( incx == 1 ) {
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
temp = x[j];
|
||||
if ( nounit ) {
|
||||
temp = temp * a[j+j*lda];
|
||||
}
|
||||
for ( i = j - 1; 0 <= i; i-- ) {
|
||||
temp = temp + a[i+j*lda] * x[i];
|
||||
}
|
||||
x[j] = temp;
|
||||
}
|
||||
} else {
|
||||
jx = kx + ( n - 1 ) * incx;
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
temp = x[jx];
|
||||
ix = jx;
|
||||
if ( nounit ) {
|
||||
temp = temp * a[j+j*lda];
|
||||
}
|
||||
for ( i = j - 1; 0 <= i; i-- ) {
|
||||
ix = ix - incx;
|
||||
temp = temp + a[i+j*lda] * x[ix];
|
||||
}
|
||||
x[jx] = temp;
|
||||
jx = jx - incx;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if ( incx == 1 ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
temp = x[j];
|
||||
if ( nounit ) {
|
||||
temp = temp * a[j+j*lda];
|
||||
}
|
||||
for ( i = j + 1; i < n; i++ ) {
|
||||
temp = temp + a[i+j*lda] * x[i];
|
||||
}
|
||||
x[j] = temp;
|
||||
}
|
||||
} else {
|
||||
jx = kx;
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
temp = x[jx];
|
||||
ix = jx;
|
||||
if ( nounit ) {
|
||||
temp = temp * a[j+j*lda];
|
||||
}
|
||||
for ( i = j + 1; i < n; i++ ) {
|
||||
ix = ix + incx;
|
||||
temp = temp + a[i+j*lda] * x[ix];
|
||||
}
|
||||
x[jx] = temp;
|
||||
jx = jx + incx;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of STRMV
|
||||
|
||||
}
|
||||
526
src/strsm.c
526
src/strsm.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <stdbool.h>
|
||||
|
||||
#include "blas.h"
|
||||
|
||||
/*
|
||||
*> \brief \b STRSM
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -177,264 +182,263 @@
|
|||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE STRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
|
||||
*
|
||||
* -- Reference BLAS level3 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
REAL ALPHA
|
||||
INTEGER LDA,LDB,M,N
|
||||
CHARACTER DIAG,SIDE,TRANSA,UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL A(LDA,*),B(LDB,*)
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MAX
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
REAL TEMP
|
||||
INTEGER I,INFO,J,K,NROWA
|
||||
LOGICAL LSIDE,NOUNIT,UPPER
|
||||
* ..
|
||||
* .. Parameters ..
|
||||
REAL ONE,ZERO
|
||||
PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
|
||||
* ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
LSIDE = LSAME(SIDE,'L')
|
||||
IF (LSIDE) THEN
|
||||
NROWA = M
|
||||
ELSE
|
||||
NROWA = N
|
||||
END IF
|
||||
NOUNIT = LSAME(DIAG,'N')
|
||||
UPPER = LSAME(UPLO,'U')
|
||||
*
|
||||
INFO = 0
|
||||
IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
|
||||
INFO = 1
|
||||
ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
|
||||
INFO = 2
|
||||
ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
|
||||
+ (.NOT.LSAME(TRANSA,'T')) .AND.
|
||||
+ (.NOT.LSAME(TRANSA,'C'))) THEN
|
||||
INFO = 3
|
||||
ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
|
||||
INFO = 4
|
||||
ELSE IF (M.LT.0) THEN
|
||||
INFO = 5
|
||||
ELSE IF (N.LT.0) THEN
|
||||
INFO = 6
|
||||
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
|
||||
INFO = 9
|
||||
ELSE IF (LDB.LT.MAX(1,M)) THEN
|
||||
INFO = 11
|
||||
END IF
|
||||
IF (INFO.NE.0) THEN
|
||||
CALL XERBLA('STRSM ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF (M.EQ.0 .OR. N.EQ.0) RETURN
|
||||
*
|
||||
* And when alpha.eq.zero.
|
||||
*
|
||||
IF (ALPHA.EQ.ZERO) THEN
|
||||
DO 20 J = 1,N
|
||||
DO 10 I = 1,M
|
||||
B(I,J) = ZERO
|
||||
10 CONTINUE
|
||||
20 CONTINUE
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Start the operations.
|
||||
*
|
||||
IF (LSIDE) THEN
|
||||
IF (LSAME(TRANSA,'N')) THEN
|
||||
*
|
||||
* Form B := alpha*inv( A )*B.
|
||||
*
|
||||
IF (UPPER) THEN
|
||||
DO 60 J = 1,N
|
||||
IF (ALPHA.NE.ONE) THEN
|
||||
DO 30 I = 1,M
|
||||
B(I,J) = ALPHA*B(I,J)
|
||||
30 CONTINUE
|
||||
END IF
|
||||
DO 50 K = M,1,-1
|
||||
IF (B(K,J).NE.ZERO) THEN
|
||||
IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
|
||||
DO 40 I = 1,K - 1
|
||||
B(I,J) = B(I,J) - B(K,J)*A(I,K)
|
||||
40 CONTINUE
|
||||
END IF
|
||||
50 CONTINUE
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
DO 100 J = 1,N
|
||||
IF (ALPHA.NE.ONE) THEN
|
||||
DO 70 I = 1,M
|
||||
B(I,J) = ALPHA*B(I,J)
|
||||
70 CONTINUE
|
||||
END IF
|
||||
DO 90 K = 1,M
|
||||
IF (B(K,J).NE.ZERO) THEN
|
||||
IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
|
||||
DO 80 I = K + 1,M
|
||||
B(I,J) = B(I,J) - B(K,J)*A(I,K)
|
||||
80 CONTINUE
|
||||
END IF
|
||||
90 CONTINUE
|
||||
100 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
*
|
||||
* Form B := alpha*inv( A**T )*B.
|
||||
*
|
||||
IF (UPPER) THEN
|
||||
DO 130 J = 1,N
|
||||
DO 120 I = 1,M
|
||||
TEMP = ALPHA*B(I,J)
|
||||
DO 110 K = 1,I - 1
|
||||
TEMP = TEMP - A(K,I)*B(K,J)
|
||||
110 CONTINUE
|
||||
IF (NOUNIT) TEMP = TEMP/A(I,I)
|
||||
B(I,J) = TEMP
|
||||
120 CONTINUE
|
||||
130 CONTINUE
|
||||
ELSE
|
||||
DO 160 J = 1,N
|
||||
DO 150 I = M,1,-1
|
||||
TEMP = ALPHA*B(I,J)
|
||||
DO 140 K = I + 1,M
|
||||
TEMP = TEMP - A(K,I)*B(K,J)
|
||||
140 CONTINUE
|
||||
IF (NOUNIT) TEMP = TEMP/A(I,I)
|
||||
B(I,J) = TEMP
|
||||
150 CONTINUE
|
||||
160 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
ELSE
|
||||
IF (LSAME(TRANSA,'N')) THEN
|
||||
*
|
||||
* Form B := alpha*B*inv( A ).
|
||||
*
|
||||
IF (UPPER) THEN
|
||||
DO 210 J = 1,N
|
||||
IF (ALPHA.NE.ONE) THEN
|
||||
DO 170 I = 1,M
|
||||
B(I,J) = ALPHA*B(I,J)
|
||||
170 CONTINUE
|
||||
END IF
|
||||
DO 190 K = 1,J - 1
|
||||
IF (A(K,J).NE.ZERO) THEN
|
||||
DO 180 I = 1,M
|
||||
B(I,J) = B(I,J) - A(K,J)*B(I,K)
|
||||
180 CONTINUE
|
||||
END IF
|
||||
190 CONTINUE
|
||||
IF (NOUNIT) THEN
|
||||
TEMP = ONE/A(J,J)
|
||||
DO 200 I = 1,M
|
||||
B(I,J) = TEMP*B(I,J)
|
||||
200 CONTINUE
|
||||
END IF
|
||||
210 CONTINUE
|
||||
ELSE
|
||||
DO 260 J = N,1,-1
|
||||
IF (ALPHA.NE.ONE) THEN
|
||||
DO 220 I = 1,M
|
||||
B(I,J) = ALPHA*B(I,J)
|
||||
220 CONTINUE
|
||||
END IF
|
||||
DO 240 K = J + 1,N
|
||||
IF (A(K,J).NE.ZERO) THEN
|
||||
DO 230 I = 1,M
|
||||
B(I,J) = B(I,J) - A(K,J)*B(I,K)
|
||||
230 CONTINUE
|
||||
END IF
|
||||
240 CONTINUE
|
||||
IF (NOUNIT) THEN
|
||||
TEMP = ONE/A(J,J)
|
||||
DO 250 I = 1,M
|
||||
B(I,J) = TEMP*B(I,J)
|
||||
250 CONTINUE
|
||||
END IF
|
||||
260 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
*
|
||||
* Form B := alpha*B*inv( A**T ).
|
||||
*
|
||||
IF (UPPER) THEN
|
||||
DO 310 K = N,1,-1
|
||||
IF (NOUNIT) THEN
|
||||
TEMP = ONE/A(K,K)
|
||||
DO 270 I = 1,M
|
||||
B(I,K) = TEMP*B(I,K)
|
||||
270 CONTINUE
|
||||
END IF
|
||||
DO 290 J = 1,K - 1
|
||||
IF (A(J,K).NE.ZERO) THEN
|
||||
TEMP = A(J,K)
|
||||
DO 280 I = 1,M
|
||||
B(I,J) = B(I,J) - TEMP*B(I,K)
|
||||
280 CONTINUE
|
||||
END IF
|
||||
290 CONTINUE
|
||||
IF (ALPHA.NE.ONE) THEN
|
||||
DO 300 I = 1,M
|
||||
B(I,K) = ALPHA*B(I,K)
|
||||
300 CONTINUE
|
||||
END IF
|
||||
310 CONTINUE
|
||||
ELSE
|
||||
DO 360 K = 1,N
|
||||
IF (NOUNIT) THEN
|
||||
TEMP = ONE/A(K,K)
|
||||
DO 320 I = 1,M
|
||||
B(I,K) = TEMP*B(I,K)
|
||||
320 CONTINUE
|
||||
END IF
|
||||
DO 340 J = K + 1,N
|
||||
IF (A(J,K).NE.ZERO) THEN
|
||||
TEMP = A(J,K)
|
||||
DO 330 I = 1,M
|
||||
B(I,J) = B(I,J) - TEMP*B(I,K)
|
||||
330 CONTINUE
|
||||
END IF
|
||||
340 CONTINUE
|
||||
IF (ALPHA.NE.ONE) THEN
|
||||
DO 350 I = 1,M
|
||||
B(I,K) = ALPHA*B(I,K)
|
||||
350 CONTINUE
|
||||
END IF
|
||||
360 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of STRSM
|
||||
*
|
||||
END
|
||||
*/
|
||||
void strsm ( char side, char uplo, char transa, char diag, int m, int n,
|
||||
float alpha, float a[], int lda, float b[], int ldb ) {
|
||||
|
||||
// -- Reference BLAS level3 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
// =====================================================================
|
||||
|
||||
// .. Local Scalars ..
|
||||
float temp;
|
||||
int i, info, j, k, nrowa;
|
||||
bool lside, nounit, upper;
|
||||
|
||||
// Test the input parameters.
|
||||
|
||||
lside = lsame ( side, 'L' );
|
||||
|
||||
if ( lside ) {
|
||||
nrowa = m;
|
||||
}
|
||||
else {
|
||||
nrowa = n;
|
||||
}
|
||||
|
||||
nounit = lsame ( diag, 'N' );
|
||||
upper = lsame ( uplo, 'U' );
|
||||
|
||||
info = 0;
|
||||
|
||||
if ( ( ! lside ) && ( ! lsame ( side, 'R' ) ) ) {
|
||||
info = 1;
|
||||
}
|
||||
else if ( ( ! upper ) && ( ! lsame ( uplo, 'L' ) ) ) {
|
||||
info = 2;
|
||||
}
|
||||
else if ( ( ! lsame ( transa, 'N' ) ) &&
|
||||
( ! lsame ( transa, 'T' ) ) &&
|
||||
( ! lsame ( transa, 'C' ) ) ) {
|
||||
info = 3;
|
||||
}
|
||||
else if ( ( ! lsame ( diag, 'U' ) ) && ( ! lsame ( diag, 'N' ) ) ) {
|
||||
info = 4;
|
||||
}
|
||||
else if ( m < 0 ) {
|
||||
info = 5;
|
||||
}
|
||||
else if ( n < 0 ) {
|
||||
info = 6;
|
||||
}
|
||||
else if ( lda < i4_max ( 1, nrowa ) ) {
|
||||
info = 9;
|
||||
}
|
||||
else if ( ldb < i4_max ( 1, m ) ) {
|
||||
info = 11;
|
||||
}
|
||||
|
||||
if ( info != 0 ) {
|
||||
xerbla ( "STRSM", info );
|
||||
return;
|
||||
}
|
||||
|
||||
// Quick return if possible.
|
||||
|
||||
if ( n == 0 ) {
|
||||
return;
|
||||
}
|
||||
|
||||
// and when alpha is 0.0.
|
||||
|
||||
if ( alpha == 0.0 ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = 0.0;
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
// Start the operations.
|
||||
|
||||
if ( lside ) {
|
||||
|
||||
// Form B := alpha*inv( a )*B.
|
||||
|
||||
if ( lsame ( transa, 'N' ) ) {
|
||||
if ( upper ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( alpha != 1.0 ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = alpha * b[i+j*ldb];
|
||||
}
|
||||
}
|
||||
for ( k = m - 1; 0 <= k; k-- ) {
|
||||
if ( b[k+j*ldb] != 0.0 ) {
|
||||
if ( nounit ) {
|
||||
b[k+j*ldb] = b[k+j*ldb] / a[k+k*lda];
|
||||
}
|
||||
for ( i = 0; i < k; i++ ) {
|
||||
b[i+j*ldb] = b[i+j*ldb] - b[k+j*ldb] * a[i+k*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( alpha != 1.0 ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = alpha * b[i+j*ldb];
|
||||
}
|
||||
}
|
||||
for ( k = 0; k < m; k++ ) {
|
||||
if ( b[k+j*ldb] != 0.0 ) {
|
||||
if ( nounit ) {
|
||||
b[k+j*ldb] = b[k+j*ldb] / a[k+k*lda];
|
||||
}
|
||||
for ( i = k + 1; i < m; i++ ) {
|
||||
b[i+j*ldb] = b[i+j*ldb] - b[k+j*ldb] * a[i+k*lda];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// Form B := alpha*inv( A' )*B.
|
||||
if ( upper ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
temp = alpha * b[i+j*ldb];
|
||||
for ( k = 0; k < i; k++ ) {
|
||||
temp = temp - a[k+i*lda] * b[k+j*ldb];
|
||||
}
|
||||
if ( nounit ) {
|
||||
temp = temp / a[i+i*lda];
|
||||
}
|
||||
b[i+j*ldb] = temp;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
for ( i = m - 1; 0 <= i; i-- ) {
|
||||
temp = alpha * b[i+j*ldb];
|
||||
for ( k = i + 1; k < m; k++ ) {
|
||||
temp = temp - a[k+i*lda] * b[k+j*ldb];
|
||||
}
|
||||
if ( nounit ) {
|
||||
temp = temp / a[i+i*lda];
|
||||
}
|
||||
b[i+j*ldb] = temp;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// Form B := alpha*B*inv( A ).
|
||||
if ( lsame ( transa, 'N' ) ) {
|
||||
if ( upper ) {
|
||||
for ( j = 0; j < n; j++ ) {
|
||||
if ( alpha != 1.0 ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = alpha * b[i+j*ldb];
|
||||
}
|
||||
}
|
||||
for ( k = 0; k < j; k++ ) {
|
||||
if ( a[k+j*lda] != 0.0 ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = b[i+j*ldb] - a[k+j*lda] * b[i+k*ldb];
|
||||
}
|
||||
}
|
||||
}
|
||||
if ( nounit ) {
|
||||
temp = 1.0 / a[j+j*lda];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = temp * b[i+j*ldb];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
for ( j = n - 1; 0 <= j; j-- ) {
|
||||
if ( alpha != 1.0 ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = alpha * b[i+j*ldb];
|
||||
}
|
||||
}
|
||||
for ( k = j + 1; k < n; k++ ) {
|
||||
if ( a[k+j*lda] != 0.0 ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = b[i+j*ldb] - a[k+j*lda] * b[i+k*ldb];
|
||||
}
|
||||
}
|
||||
}
|
||||
if ( nounit ) {
|
||||
temp = 1.0 / a[j+j*lda];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = temp * b[i+j*ldb];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Form B := alpha*B*inv( A' ).
|
||||
|
||||
else {
|
||||
if ( upper ) {
|
||||
for ( k = n - 1; 0 <= k; k-- ) {
|
||||
if ( nounit ) {
|
||||
temp = 1.0 / a[k+k*lda];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+k*ldb] = temp * b[i+k*ldb];
|
||||
}
|
||||
}
|
||||
for ( j = 0; j < k; j++ ) {
|
||||
if ( a[j+k*lda] != 0.0 ) {
|
||||
temp = a[j+k*lda];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = b[i+j*ldb] - temp * b[i+k*ldb];
|
||||
}
|
||||
}
|
||||
}
|
||||
if ( alpha != 1.0 ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+k*ldb] = alpha * b[i+k*ldb];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else {
|
||||
for ( k = 0; k < n; k++ ) {
|
||||
if ( nounit ) {
|
||||
temp = 1.0 / a[k+k*lda];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+k*ldb] = temp * b[i+k*ldb];
|
||||
}
|
||||
}
|
||||
for ( j = k + 1; j < n; j++ ) {
|
||||
if ( a[j+k*lda] != 0.0 ) {
|
||||
temp = a[j+k*lda];
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+j*ldb] = b[i+j*ldb] - temp * b[i+k*ldb];
|
||||
}
|
||||
}
|
||||
}
|
||||
if ( alpha != 1.0 ) {
|
||||
for ( i = 0; i < m; i++ ) {
|
||||
b[i+k*ldb] = alpha * b[i+k*ldb];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
// End of STRSM
|
||||
|
||||
}
|
||||
47
src/xerbla.c
47
src/xerbla.c
|
|
@ -1,3 +1,8 @@
|
|||
#include <stdio.h>
|
||||
|
||||
#include "errquit.h"
|
||||
|
||||
/*
|
||||
*> \brief \b XERBLA
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
|
|
@ -56,31 +61,17 @@
|
|||
*> \ingroup aux_blas
|
||||
*
|
||||
* =====================================================================
|
||||
SUBROUTINE XERBLA( SRNAME, INFO )
|
||||
*
|
||||
* -- Reference BLAS level1 routine --
|
||||
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
CHARACTER*(*) SRNAME
|
||||
INTEGER INFO
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC LEN_TRIM
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
WRITE( *, FMT = 9999 )SRNAME( 1:LEN_TRIM( SRNAME ) ), INFO
|
||||
*
|
||||
STOP
|
||||
*
|
||||
9999 FORMAT( ' ** On entry to ', A, ' parameter number ', I2, ' had ',
|
||||
$ 'an illegal value' )
|
||||
*
|
||||
* End of XERBLA
|
||||
*
|
||||
END
|
||||
*/
|
||||
void xerbla( char *srname, int info ) {
|
||||
|
||||
// -- Reference BLAS level1 routine --
|
||||
// -- Reference BLAS is a software package provided by Univ. of Tennessee, --
|
||||
// -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
|
||||
printf("xerbla:double: lapack error\n");
|
||||
printf(" ** On entry to %s parameter number %d had an illegal value", srname, info);
|
||||
exit(1);
|
||||
|
||||
// End of XERBLA
|
||||
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue