HP Forums
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)<SIZE(Y)) 
 BEGIN
  RET (VECX);
 END;

 //CAXPY - y = a*x + y
 EXPORT CAXPY(SCL,VECX,VECY)    //   SUBROUTINE CAXPY(N,CA,CX,INCX,CY,INCY)
 //*>    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)<SIZE(Y)) 
 BEGIN
  RET (VECX);//IDENTICAL
 END;

 //SAXPY - y = a*x + y //const*VEC+VEC
 EXPORT SAXPY(SCL,VECX,VECY) // SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)
 BEGIN
  RETURN CAXPY(SCL,VECX,VECY);//IDENTICAL (ASSUMING REAL VECS)
  //RETURN {(SCL*VECX+VECY)};//IDENTICAL
 END;

 //SDOT - dot product //of 2 vectors
 EXPORT SDOT(VECX,VECY) // REAL SDOT(N,SX,INCX,SY,INCY)
 BEGIN
  //PRINT({VECX,VECY});
  RETURN CDOTU(VECX,VECY);//IDENTICAL (ASSUMING VECS ARE REAL)
  //RETURN dot(VECX,VECY);//IDENTICAL
 END;

 //SDSDOT - dot product with extended precision accumulation
 EXPORT SDSDOT(SB,VECX,VECY) //REAL FUNCTION SDSDOT(N,SB,SX,INCX,SY,INCY)
 //THE ORIGINAL FORTRAN IMPLEMENTS DOT USING DOUBLES TO IMPROVE PRECISION
 //PPL: NO EXTENDED PRECISION (iust making parameters of dot exact doesnt help)
 BEGIN
  RETURN SB+SDOT(VECX,VECY);
 END;

 //SNRM2 - Euclidean norm
 EXPORT SNRM2(VECX) //REAL FUNCTION SNRM2(N,X,INCX)
 BEGIN
  LOCAL EUC:=CAS.l2norm(VECX); 
  RETURN IFTE(RAN,approx(EUC),EUC);
 END;
 //These currently return exact if RAN=0
 //SCNRM2- Euclidean norm of complex
 EXPORT SCNRM2(CVECX)
 BEGIN
  LOCAL EUC:=CAS.l2norm(CVECX);
  RETURN IFTE(RAN,approx(EUC),EUC);
 END;

 //SASUM - sum of absolute values
 EXPORT SASUM(VEC) //       REAL FUNCTION SASUM(N,SX,INCX)
 BEGIN 
  //RETURN ΣLIST(ABS(mat2list(VEC)));
  RETURN CAS.l1norm(VEC);
  //RETURN CSASUM(VEC);//IDENTICAL
 END;

 //ISAMAX - index of max abs value
 EXPORT ISAMAX(VEC) //INTEGER FUNCTION ISAMAX(N,SX,INCX)
 //INDEX OF MAX ABS REAL
 //INDEX OF MAX ABS (REALPART+IMPART) Z
 //NB THE CODE IS OPTIMISED ASSUMING REAL.
 //IT WILL GIVE DIFFERENT RESULTS THAN ICAMAX IF GIVEN COMPLEX.
 //IF YOU WANT IDENTICAL RESULTS CALL ICAMAX
 BEGIN 
  LOCAL STRIDE:=1;
  LOCAL LVC:=mat2list(VEC);
  LOCAL LBOTH:=ABS(LVC);//OPTIMISED FOR REAL.
  //ABS(RE(LVC))+ABS(IM(LVC));//COMPLEX

  //IF STRIDE<0 THEN NONSTANDARD
  // RETURN ICAMAX(list2mat(REVERSE(LVC)),−STRIDE);//TBD
  //END;
  IF STRIDE>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.