BLAS-Prime: Basic Linear Algebra Subprograms - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP Prime Software Library (/forum-15.html) +--- Thread: BLAS-Prime: Basic Linear Algebra Subprograms (/thread-13613.html) BLAS-Prime: Basic Linear Algebra Subprograms - StephenG1CMZ - 09-08-2019 11:11 PM I have been comparing the functionality of the common Basic Linear Algebra Subprograms and those operations built into the Prime. In doing so, I have implemented some BLAS functions in PPL, based on the Fortran reference code that is unoptimised. Or using a built-in if obvious. RE: BLAS-Prime: Basic Linear Algebra Subprograms - StephenG1CMZ - 09-08-2019 11:13 PM V 0.1 implements most of BLAS Level 1 functions in single precision. Main code: Code: ```   //CUSTOMISE  LOCAL RABL:=1; //RETURN AS BLAS OR LIST   //TBD:BETTER EXAMPLE   //SET RAL TO 1: ENABLES THE BLAS-LIKE SYNTAX     //sto(SSWAP(XX,YY),{'XX','YY'}); //XX AND YY UPDATE      //a list of parameters can be updated   //SET RAL:=0: USE THE PPL SYNTAX     //RESULTS:=SSWAP(XX,YY);     //XX:=RESULTS(1);YY:=RESULTS(2);...ETC   LOCAL RAN:=1; //RESULT APPROXIMATE NUMERIC    //1-(RESULT APPROX NUMERIC: DEFAULT)    //0-(RESULT PPL EXACT WHEN AVAILABLE FROM IMPLEMENT)   //END CUSTOM  EXPORT ABOUT()  BEGIN   PRINT();   PRINT("BLAS-Prime V0.1 © StephenG1CMZ");   PRINT("Based on BLAS Reference Fortran code (not optimised)");   PRINT("http://www.netlib.org/blas/#_reference_blas_version_3_8_0");   PRINT("This PPL implementation is not verified ");   PRINT("V0.1: mimics BLAS Level 1");   PRINT("No N or stride");   PRINT("Some Givens rotations TBD");  END;  //PPL SUPPORT ROUTINES  RET(PARAM)  BEGIN   RETURN IFTE(RABL,{PARAM},PARAM);  END; EXPORT SFPINFO () //ENVIRONMENT BEGIN END;  EXPORT MYCONJ(VEC)  BEGIN   LOCAL CNJVEC;   //CNJVC:=list2mat(conj(mat2list(VEC)));   LOCAL LVEC:=mat2list(VEC);    LOCAL LCNJVEC:=(RE(LVEC)-IM(LVEC)*);//CONJ NEGATES IM   CNJVEC:=list2mat(LCNJVEC,SIZE(LCNJVEC));   //PRINT("CNJ: "+CNJVEC);   RETURN CNJVEC(1);//AN ARTIFACT OF SYNTAX ON VECTORS  END;  //BLAS  EXPORT BLAS_L1_vector() //HEADING  BEGIN  END;  //BLAS LEVEL 1 VECTOR OPERATIONS  //Complex Single    //CROTG - setup Givens rotation  EXPORT CSROTGTBD()  BEGIN  END;  //CSROT - apply Givens rotation  EXPORT CSROTTBD()  BEGIN  END;  //CSWAP - swap x and y  EXPORT CSWAP(VECX,VECY)  //call: sto(SSWAP(XX,YY),{'XX','YY'});  BEGIN   RETURN {VECY,VECX};//IDENTICAL  END;  //CSCAL - x = a*x  EXPORT CSCAL(SCL,VECX)    //   SUBROUTINE CSCAL(N,CA,CX,INCX)  //*>    CSCAL scales a vector by a constant.  //*> \param[in] N  //*> \param[in] CA       CA is COMPLEX  //*> \param[in,out] CX  CX Is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )  BEGIN   RET  (SCL*VECX);//IDENTICAL   END;  //CSSCAL - x = a*x  EXPORT CSSCAL(SCL,VECX)     //  SUBROUTINE CSSCAL(N,SA,CX,INCX)  //*>    CSSCAL scales a complex vector by a real constant.  //*> \param[in] N  //*> \param[in] SA         SA is REAL  //*> \param[in,out] CX     CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )  BEGIN   RET (SCL*VECX);//IDENTICAL  END;    //CCOPY - copy x into y  EXPORT CCOPY(VECX,VECY)   //    SUBROUTINE CCOPY(N,CX,INCX,CY,INCY)  //*>    CCOPY copies a vector x to a vector y.  //*> \param[in] N  //*> \param[in] CX    CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )  //*> \param[out] CY   CY is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCY ) )  //NOTE:MY IMPLEMENT CURRENTLY REPLACES Y BY X (LEAVING IT SHORTER WHEN SIZE(X)    CAXPY constant times a vector plus a vector.  //*> \param[in] N  //*> \param[in] CA    CA is COMPLEX  //*>           On entry, CA specifies the scalar alpha.  //*> \param[in] CX    CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )  //*> \param[in,out] CY   CY is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCY ) )  BEGIN   RETURN {(SCL*VECX+VECY)};//IDENTICAL  END;  //CDOTU - dot product  EXPORT CDOTU(VECX,VECY)    //   COMPLEX FUNCTION CDOTU(N,CX,INCX,CY,INCY)  //      COMPLEX CX(*),CY(*)  //*> CDOTU forms the dot product of two complex vectors  //*>      CDOTU = X^T * Y  //*> \param[in] CX    CX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )  //*> \param[in] CY    CY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )  BEGIN   //LOCAL II;   //LOCAL CTMP:=0;   //LOCAL DIMX:=DIM(VECX);   //LOCAL DIMY:=DIM(VECY);   //LOCAL NNN:=MIN({DIMX(1),DIMX(1)});   //MSGBOX(DIM(VECX));   //FOR II:=1 TO NNN DO   // CTMP:=CTMP+VECX(II)*VECY(II);   //END;   //RETURN CTMP;      RETURN dot(VECX,VECY);//IDENTICAL  END;  //CDOTC - dot product, conjugating the first vector  EXPORT CDOTC(VECX,VECY)   //    COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY)  //*       COMPLEX CX(*),CY(*)  //*> CDOTC forms the dot product of two complex vectors  //*>      CDOTC = X^H * Y  //*> \param[in] CX        CX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX ) )  //*> \param[in] CY        CY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )  BEGIN   //PRINT({VECX,VECY});   //PRINT({MYCONJ(VECX),VECY});   RETURN CDOTU(MYCONJ(VECX),VECY);  END;  //SCASUM - sum of absolute values  EXPORT SCASUM(VEC) //  REAL FUNCTION SCASUM(N,CX,INCX)  //*>    SCASUM takes the sum of the (|Re(.)| + |Im(.)|)'s of a complex vector and  //*>    returns a single precision result.  //*> \param[in,out] CX       CX is COMPLEX array, dimension ( 1 + ( N - 1 )*abs( INCX ) )  //Query: CX is inout but unchanged  //sum(abs) re and im separately  BEGIN   //RETURN SASUM(RE(mat2list(VEC)))+SASUM(IM(mat2list(VEC)))*;   RETURN CAS.l1norm(VEC);   //RETURN SASUM(VEC);//IDENTICAL  END;  //ICAMAX - index of max abs value  EXPORT ICAMAX(VC)  BEGIN  LOCAL STRIDE:=1;  //ICAMAX and IZAMAX find the first element xk, where k is defined as the smallest index k, such that:  //|ak|+|bk| = max{|aj|+|bj| for j = 1, n}  //where xk = (ak, bk)  //By specifying a positive or negative stride for vector x, the first or last occurrence, respectively, is found in the array. The position i, returned as the value of the function, is always figured relative to the location specified in the calling sequence for vector x (in argument x). Therefore, depending on the stride specified for incx, i has the following values:  //For incx ≥ 0, i = k  //For incx < 0, i = n-k+1   LOCAL LVC:=mat2list(VC);   LOCAL LBOTH:=ABS(RE(LVC))+ABS(IM(LVC));   //IF STRIDE<0 THEN NONSTANDARD   // RETURN ICAMAX(list2mat(REVERSE(LVC)),−STRIDE);//TBD   //END;   IF STRIDE>1 THEN    //TBD;   END;   RETURN POS(LBOTH,MAX(LBOTH));  END;  //Real Single  //SROTG - setup Givens rotation  EXPORT SROTGTBD(SA,SB) // SUBROUTINE SROTG(SA,SB,C,S)*     .. Scalar Arguments ..  //REAL C,S,SA,SB  BEGIN   LOCAL CC,SS;   MSGBOX("TBD");   RETURN {CC,SS};  END;  //SROTMG - setup modified Givens rotation  EXPORT SROTMGTBD(SD1,SD2,SX1,SY1) //       SUBROUTINE SROTMG(SD1,SD2,SX1,SY1,SPARAM)  //*> \param[in,out] SD1  SD1 is REAL  //*> \param[in,out] SD2  SD2 is REAL  //*> \param[in,out] SX1  SX1 is REAL  //*> \param[in]     SY1  SY1 is REAL  //*> \param[out] SPARAM  SPARAM is REAL array, dimension (5)   //*>     SPARAM(1)=SFLAG   //*>     SPARAM(2)=SH11   //*>     SPARAM(3)=SH21   //*>     SPARAM(4)=SH12   //*>     SPARAM(5)=SH22  BEGIN   LOCAL SD1OUT,SD2OUT,SX1OUT;   LOCAL PARAM:={};   MSGBOX("TBD");   RETURN {SD1OUT,SD2OUT,SX1OUT,PARAM};  END;  //SROT - apply Givens rotation  EXPORT SROTTBD(VECSX,VECSY,CC,SS) //SUBROUTINE SROT(N,SX,INCX,SY,INCY,C,S)  //*> \param[in,out] SX   SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX )  //*> \param[in]     INCX INCX is INTEGER  //*> \param[in,out] SY   SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )  //*> \param[in]     INCY INCY is INTEGER  //*> \param[in]     C    C is REAL  //*> \param[in]     S    S is REAL  BEGIN   LOCAL SX,SY;   MSGBOX("TBD");   RETURN {SX,SY};  END;  //SROTM - apply modified Givens rotation  EXPORT SROTMTBD(VECSX,VECSY,SPARAM) //SUBROUTINE SROTM(N,SX,INCX,SY,INCY,SPARAM)  //*> \param[in]       N   N is INTEGER  //*>         number of elements in input vector  //*> \param[in,out]   SX SX is REAL array, dimension ( 1 + ( N - 1 )*abs( INCX )  //*> \param[in,out]   SY SY is REAL array, dimension ( 1 + ( N - 1 )*abs( INCY ) )  //*> \param[in]  SPARAM  SPARAM is REAL array, dimension (5)   //*>     SPARAM(1)=SFLAG   //*>     SPARAM(2)=SH11   //*>     SPARAM(3)=SH21   //*>     SPARAM(4)=SH12   //*>     SPARAM(5)=SH22  BEGIN   LOCAL SX,SY;   RETURN {SX,SY};  END;  //SSWAP - swap x and y  EXPORT SSWAP(VECX,VECY) //   SUBROUTINE SSWAP(N,SX,INCX,SY,INCY)  //call: sto(SSWAP(XX,YY),{'XX','YY'});  BEGIN   RETURN {VECY,VECX};//IDENTICAL  END;  //SSCAL - x = a*x  EXPORT SSCAL(SCL,VEC) //SUBROUTINE SSCAL(N,SA,SX,INCX)  //call: sto(SSCAL(SCL,VEC),{'VEC'});  BEGIN   RET (SCL*VEC);//IDENTICAL  END;  //SCOPY - copy x into y  EXPORT SCOPY(VECX,VECY) // SUBROUTINE SCOPY(N,SX,INCX,SY,INCY)   //NOTE:MY IMPLEMENT CURRENTLY REPLACES Y BY X (LEAVING IT SHORTER WHEN SIZE(X)1 THEN    //TBD;   END;   RETURN POS(LBOTH,MAX(LBOTH));   //RETURN ICAMAX(VEC);//UNOPTIMISED  END;  //EXTRAS: NOT BLAS STANDARD  //REAL SINGLE  EXPORT ISMAX (VC)  //INDEX MAXIMUM (CF ISAMAX)  //OK REAL  BEGIN   LOCAL LVC:=mat2list(VC);   LOCAL MX:=MAX(LVC);   RETURN POS(LVC,MX);   END;   EXPORT ISMIN (VC)  //INDEX MINIMUM (CF ISMAX ISAMAX)  //OK REAL  BEGIN   LOCAL LVC:=mat2list(VC);   LOCAL MN:=MIN(LVC);   RETURN POS(LVC,MN);   END;   EXPORT BLAS()  BEGIN   //PPL TYPE DOES NOT DISTINGUISH VECTORS AND MATRICES   //PRINT(TYPE([1]));   //PRINT(TYPE([[1],[2]]));   END;``` Some Givens rotations: Code: ```    EXPORT CABS(NUM)  //FTN COMPLEX ABS  BEGIN   RETURN ABS(NUM);   //RETURN SQRT(RE(NUM)^2+IM(NUM)^2);  END;  EXPORT SIGNFTN(AA,BB)  //FTN COPYSIGN RETURNS SGN(BB)*AA   //IF B ZERO ABS(A)  //RETURN {1,0,−1} (SOME IMPLEMENTS ALLOW −0 RETURN)  BEGIN   RETURN IFTE(BB<0,-ABS(AA),ABS(AA));  END;  EXPORT COPYSIGN(AVALUE,ASIGN)  BEGIN   RETURN IFTE(ASIGN<0,-ABS(AVALUE),ABS(AVALUE));   //RETURN SIGNFTN(AVALUE,ASIGN); END;  EXPORT CONJG(ZZ)  //CONJUGATE OF A COMPLEX NUMBER  //FTN  BEGIN   RETURN RE(ZZ)-IM(ZZ)*;  END;  EXPORT CROTG(CXA,CXB) // *       SUBROUTINE CROTG(CA,CB,C,S)  //*       COMPLEX CA,CB,S  //*       REAL C  //CROTG determines a complex Givens rotation.  //*> \param[in] CA       CA is COMPLEX  //*> \param[in] CB       CB is COMPLEX  //*> \param[out] C        C is REAL  //*> \param[out] S        S is COMPLEX  //call: sto(CROTG(AA,BB),{'AA','BB','CC','SS'})  BEGIN   LOCAL CA:=CXA; //LOCAL CB:=CXB;//INOUTS   LOCAL CC,SS;     // COMPLEX SS     // REAL CC       //*     .. Local Scalars ..    LOCAL ALPHA; //COMPLEX     LOCAL NORM,SCALE; //REAL    //*     .. Intrinsic Functions ..    //    INTRINSIC CABS,CONJG,SQRT     //PRINT({"CROTG IN:",CA,CXB});       IF (CABS(CA)==0.) THEN          CC:= 0.;          SS:= 1.0; //(1.,0.);          CA:= CXB;       ELSE          SCALE:= CABS(CA) + CABS(CXB);          //SCALE:=1;//SIMPLER          NORM:= SCALE*SQRT((CABS(CA/SCALE))^2+ (CABS(CXB/SCALE))^2);          ALPHA:= CA/CABS(CA);          CC:= CABS(CA)/NORM;          SS:= ALPHA*CONJG(CXB)/NORM;          CA:= ALPHA*NORM;       END;//IF;       //RETURN ROUND({CA,CC,SS},2);//SIMPLER VALUES IN TESTS       RETURN {CA,CC,SS}; //RETURN 2 COMPLEXES AND CC REAL? (CB UNCHANGED)  END;  EXPORT CSROT(VECX,VECY,CS,SN) //       SUBROUTINE CSROT( N, CX, INCX, CY, INCY, CS, SN)  //*       .. Scalar Arguments ..  //*       INTEGER           INCX, INCY, N  //*       REAL              CS, SN  //*       .. Array Arguments ..  //*       COMPLEX           CX( * ), CY( * )  //*> CSROT applies a plane rotation, where the cos and sin (c and s) are real  //*> and the vectors cx and cy are complex.  //*> \param[in,out] CX     CX is COMPLEX array,  //*> \param[in,out] CY     CY is COMPLEX array,   //*> \param[in] CS         CS is REAL  //*>           On entry, C specifies the cosine, cos.  //*> \param[in] SN         SN is REAL  BEGIN   LOCAL VECXOUT,VECYOUT;      VECXOUT:=CS*VECX+SN*VECY;   VECYOUT:=CS*VECY-SN*VECX;      RETURN ROUND({VECXOUT,VECYOUT},3);//TESTS   RETURN {VECXOUT,VECYOUT};  END;  EXPORT SROTG(SA,SB) //  SUBROUTINE SROTG(SA,SB,C,S)  //  REAL C,S,SA,SB  //*>          SA is REAL  //*>          SB is REAL  //*> \param[out] C     C is REAL  //*> \param[out] S     S is REAL  //call: sto(SROTG(AA,BB),{'AA','BB','CC','SS'})  //*>    SROTG construct givens plane rotation.   //A real Givens plane rotation is constructed for values a and b by computing values for r, c, s, and z, where:  //Real Givens Plane Rotation Math Graphic    BEGIN   LOCAL RR,ROE,SCALE,ZZ;   LOCAL CC;    LOCAL SS;       ROE:=IFTE(ABS(SA)>ABS(SB),SA,SB);       SCALE:= ABS(SA) + ABS(SB);//AVOID LOP       IF (SCALE==0.0) THEN          CC:= 1.0;          SS:= 0.0;          RR:= 0.0;          ZZ:= 0.0;       ELSE          RR:=SCALE*SQRT((SA/SCALE)^2+ (SB/SCALE)^2); //AVOID LOP           //FTN:SCALING TO AVOID LOP           //PPL:EXAMPLES BETTER WITHOUT (ALT. ROUND OUTPUTS)          //RR:=SQRT((SA)^2+ (SB)^2);          RR:=COPYSIGN(1.0,ROE)*RR; //COPYSIGN           //RR:=SIGN(ROE)*RR;          CC:= SA/RR;          SS:= SB/RR;          ZZ:= 1.0;          IF (ABS(SA)>ABS(SB)) THEN ZZ:= SS; END;          IF (ABS(SB)≥ABS(SA) AND CC≠0.0) THEN ZZ:= 1.0/CC; END;       END;// IF       //SA := RR;//RR AND ZZ EXPECTED TO REPLACE SA SB       //SB := ZZ;       //RETURN ROUND({RR,ZZ,CC,SS},3);//SIMPLER VALUES IN TESTS       RETURN {RR,ZZ,CC,SS}; //RETURN 4 REALS  END;  EXPORT EXGIVENS()  BEGIN   LOCAL CS:=0.5;   LOCAL SN:=(√(3.0))/2.0;   PRINT();   PRINT({SROTG(0.0,0.0),"=",{0,0,1,0}});   PRINT({SROTG(0.0,2.0),"=",{2,1,0,1}});   PRINT({SROTG(6.0,−8.0),"=",{−10,−1.667,−0.6,0.8}});    PRINT({SROTG(8.0,6.0),"=",{10,0.6,0.8,0.6}});   PRINT({CROTG(0.0,1.0),"=", "A=(1.0, 0.0) C =  0.0 S =  (1.0, 0.0)"});    PRINT({CROTG((3.0+ 4.0*),(4.0+6.0*)),"=", "  A =  (5.26, 7.02) C  =  0.57 S  =  (0.82, -0.05) "});   PRINT({CSROT([(1.0+ 2.0*),(2.0+ 3.0*), (3.0+ 4.0*)],[(-1.0+ 5.0*),(-2.0+ 4.0*),(-3.0+ 3.0*)],CS,SN)});   PRINT({"=","X=((-0.366, 5.330), (-0.732, 4.964), (-1.098, 4.598)) Y=((-1.366, 0.768),(-2.732, -0.598),(-4.098, -1.964))"})  END;  EXPORT GIVENS()  BEGIN   //A VECTOR MUST BE <20000   //BUT SOME ROUTINES USE MATLIST SO 10000   LOCAL CS:=0.5;   LOCAL SN:=(√3)/2;   LOCAL VECX:=MAKEMAT(10*J/1000,1000);   LOCAL VECY:=MAKEMAT(10*J/1000,1000);   PRINT();   PRINT("GO "+DIM(VECX));   PRINT({TEVAL(CROTG((3.0+ 4.0*),(4.0+6.0*)))});   PRINT({TEVAL(CSROT(VECX,VECY,CS,SN))});  END;``` Vectors in PPL can have 20 000-1 elements, but some functions are implemented using mat2list, and lists only have 10 000 elements. There is currently no error checking. Note: Some functions may round results. This is simply for clarity in displaying test results.