C ALGORITHM 830, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 30, NO. 1, March, 2004, P. 86-- 94. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Fortran90/ # Fortran90/BenchMark/ # Fortran90/BenchMark/DBenchmark.f90 # Fortran90/BenchMark/Dbenchmark.res # Fortran90/BenchMark/Givens_Rotations.f90 # Fortran90/BenchMark/Makefile # Fortran90/BenchMark/SBenchmark.f90 # Fortran90/BenchMark/Sbenchmark.res # Fortran90/BenchMark/data_dp # Fortran90/BenchMark/data_sp # Fortran90/Drivers/ # Fortran90/Drivers/DXrot.f90 # Fortran90/Drivers/Makefile # Fortran90/Drivers/SXrot.f90 # Fortran90/Drivers/d_rot_battery.f90 # Fortran90/Drivers/d_rotm_battery.f90 # Fortran90/Drivers/data1 # Fortran90/Drivers/data3 # Fortran90/Drivers/ddate_res # Fortran90/Drivers/gsrot_battery.f90 # Fortran90/Drivers/gsrotm_battery.f90 # Fortran90/Drivers/res1 # Fortran90/Drivers/res10 # Fortran90/Drivers/res11 # Fortran90/Drivers/res12 # Fortran90/Drivers/res2 # Fortran90/Drivers/res3 # Fortran90/Drivers/res4 # Fortran90/Drivers/res5 # Fortran90/Drivers/res6 # Fortran90/Drivers/res7 # Fortran90/Drivers/res8 # Fortran90/Drivers/res9 # Fortran90/Drivers/rot_battery.f90 # Fortran90/Drivers/rotm_battery.f90 # Fortran90/Drivers/s_rot_battery.f90 # Fortran90/Drivers/s_rotm_battery.f90 # Fortran90/Drivers/sdate_res # Fortran90/Src/ # Fortran90/Src/Givens_Rotations.f90 # Fortran90/Src/Import_Kinds.f90 # This archive created: Tue Apr 6 14:25:45 2004 export PATH; PATH=/bin:$PATH if test ! -d 'Fortran90' then mkdir 'Fortran90' fi cd 'Fortran90' if test ! -d 'BenchMark' then mkdir 'BenchMark' fi cd 'BenchMark' if test -f 'DBenchmark.f90' then echo shar: will not over-write existing file "'DBenchmark.f90'" else cat << "SHAR_EOF" > 'DBenchmark.f90' FUNCTION elaptime(old_time) ! Use this timing function in pairs: ! TS = ELAPTIME(TICKS); { Time something}; TS=ELAPTIME(TICKS) ! The first call starts the clock by recording TICKS. ! (Ignore TS the first call.) ! The second call updates TICKS and computes TS, in Seconds. ! This routine allows for 1-clock roll over. ! This is a Fortran 90 standard routine. It may return 0. if the ! underlying processor has no clock. ! .. ! .. Function Return Value .. REAL :: elaptime ! .. ! .. Scalar Arguments .. REAL, INTENT (INOUT) :: old_time ! .. ! .. Local Scalars .. REAL :: newtime ! .. ! .. Intrinsic Functions .. INTRINSIC cpu_time ! .. CALL cpu_time(newtime) elaptime = newtime-old_time old_time=newtime END FUNCTION elaptime PROGRAM main_timer ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: zero = 0.0_wp LOGICAL, PARAMETER :: details = .FALSE. INTEGER, PARAMETER :: mtests=4 character (3), parameter :: month(12) = (/'Jan', 'Feb',& 'Mar','Apr','May','Jun','Jul','Aug','Sep','Oct', & 'Nov', 'Dec'/) ! .. ! .. Local Scalars .. INTEGER :: i, ierror, m, mon, n, ntimes, repeat_timings LOGICAL :: success CHARACTER (64) :: identification, file_name character (8):: date character (11):: time ! .. ! .. Local Arrays .. REAL (wp) :: maxtim(mtests), mintim(mtests) REAL (wp) :: mtime(mtests), sd(mtests) REAL (wp), ALLOCATABLE :: times(:,:) LOGICAL :: passed(mtests) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: test_pass ! .. ! .. External Subroutines .. EXTERNAL givens_test, output_run_details ! .. ! .. Intrinsic Functions .. INTRINSIC minloc, minval, real, sum, trim ! .. ! Setting details to true will generate detailed timings of ! all methods for each run. WRITE (*,fmt='(/1x,A)') & ' Benchmark Givens Transformations '// & &'on 2N by N random matrices:' WRITE (*,fmt='(1x,A)',advance='NO') ' ENTER the number of COLUMNS: >' READ (*,*,err=10,end=7) n WRITE (*,fmt='(1x,A)',advance='NO') & ' ENTER the number of times to test: >' READ (*,*,err=10) ntimes WRITE (*,fmt='(1x,A)',advance='NO') & ' ENTER platform identification comment: >' READ (*,fmt='(A64)',err=10) identification WRITE (*,fmt='(1x,A)',advance='NO') & ' ENTER how many times to repeat the whole test: >' READ (*,*) repeat_timings WRITE (*,fmt='(1x,A)',advance='NO') & ' File name for results: >' READ (*,'(a)')file_name open(unit=4, file=file_name, status='unknown', iostat=ierror) IF (ierror/=0)THEN WRITE (*,'(" Error attempting to open results file: ",a)') & & file_name STOP END IF ! ! Print header information call date_and_time(date,time) read(date(5:6),'(i2)')mon WRITE (4,fmt='(15x,''Givens Transformation Benchmark Program'')') WRITE (4,fmt='(22x,''Double Precision Version''/)') WRITE (4,fmt='(15x,''Run started: '',a2,'':'',a2,'' on '', & &a2,'' '',a3,'' '',a4)') time(1:2), time(3:4), & &date(7:8), month(mon), date(1:4) m = 2*n ALLOCATE (times(mtests,repeat_timings),STAT=ierror) IF (ierror>0) THEN WRITE (4,'(" Error attempting to allocate timing array")') STOP END IF success = .TRUE. DO i = 1, repeat_timings CALL givens_test(n,ntimes,times(1:mtests,i),passed) success = success .AND. test_pass(passed) IF (details) THEN CALL output_run_details(times(1:mtests,i),success) END IF END DO ! Compute average time mtime = zero DO i = 1, repeat_timings mtime = mtime + times(:,i) END DO mtime = mtime/real(repeat_timings,wp) ! and standard deviations sd = zero DO i = 1, repeat_timings sd = sd + (times(:,i)-mtime)**2 END DO DO i = 1, mtests maxtim(i) = maxval(times(i,:)) mintim(i) = minval(times(i,:)) ENDDO ! Output synopsis of runs i = sum(minloc(mtime)) WRITE (4,'(15x,"Platform id: ",a/)') trim(identification) WRITE (4,'(15x,"Least squares problem size: ",i4," by ",i4)') m, n WRITE (4,'(15x,"Problem solved ",i5," times per run")')ntimes WRITE (4,'(15x,"Timings averaged over ",i3," complete runs"/)')& &repeat_timings WRITE (4,fmt='(25x,"Timing Details")') WRITE (4,fmt='(25x,"--------------"/)') WRITE (4,fmt='(28x,"New MG",4x,"New SG",3x,"Orig MG",3x,"Orig SG")') 100 format(a22,":",8f10.2) WRITE (4,100)"Average Times", mtime WRITE (4,100)"Standard Deviations", sd WRITE (4,100)"Normalized", mtime/minval(mtime) write (4,'()') WRITE (4,100)"Max Times", maxtim WRITE (4,100)"Min Times", mintim WRITE (4,100)"Normalized Min Times", mintim/minval(mintim) WRITE (4,fmt='(/5x,A)',advance='NO') & & " *** Optimal timings obtained using " SELECT CASE (i) CASE (1) WRITE (4,'("new MG")') CASE (2) WRITE (4,'("new SG")') CASE (3) WRITE (4,'("original Blas-1 MG")') CASE (4) WRITE (4,'("original Blas-1 SG")') END SELECT WRITE (4,'(/)') DEALLOCATE (times) 7 STOP 10 CONTINUE WRITE (4,'(" Error reading input")') END PROGRAM main_timer SUBROUTINE compute_sol(x,a) ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Array Arguments .. REAL (wp), INTENT (IN) :: a(:,:) REAL (wp), INTENT (INOUT) :: x(:) ! .. ! .. Local Scalars .. INTEGER :: i, j, n ! .. ! .. Intrinsic Functions .. INTRINSIC size ! .. n = size(x) DO j = n, 1, -1 x(j) = x(j)/a(j,j) DO i = 1, j - 1 x(i) = x(i) - a(i,j)*x(j) END DO END DO END SUBROUTINE compute_sol FUNCTION test_pass(passed) RESULT (success) ! .. Function Return Value .. LOGICAL :: success ! .. ! .. Parameters .. INTEGER, PARAMETER :: maxval = 4 ! .. ! .. Array Arguments .. LOGICAL, INTENT (IN) :: passed(maxval) ! .. ! .. Local Scalars .. INTEGER :: i LOGICAL :: ltemp ! .. success = .TRUE. DO i = 1, maxval, 2 ltemp = passed(i) .AND. passed(i+1) success = success .AND. ltemp IF ( .NOT. ltemp) THEN SELECT CASE ((i+1)/2) CASE (1) WRITE (4,*) 'Error in test: new MG then SG' CASE (2) WRITE (4,*) 'Error in test: original Blas-1-1 MG then SG' CASE DEFAULT WRITE (4,*) 'Illegal value of i ', i, ' in test_pass' END SELECT END IF END DO END FUNCTION test_pass SUBROUTINE output_run_details(times,success) ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Parameters .. INTEGER, PARAMETER :: maxval = 4 ! .. ! .. Scalar Arguments .. LOGICAL :: success ! .. ! .. Array Arguments .. REAL (wp) :: times(maxval) ! .. ! .. Local Scalars .. INTEGER :: i ! .. IF (success) THEN WRITE (4,'(" All tests successfully ran")') ELSE WRITE (4,'(" IGNORE TIMINGS: One or more tests failed")') END IF DO i = 1, maxval, 2 SELECT CASE ((i+1)/2) CASE (1) WRITE (4,'(1x,A,1x,2F7.2)') 'Timings: new MG then SG:' & , times(i), times(i+1) CASE (2) WRITE (4,'(1x,A,1x,2F7.2)') & &'Timings: original Blas-1-1 MG then SG:',& ×(i), times(i+1) CASE DEFAULT WRITE (4,*) 'Illegal value of i ', i, ' in test_pass' END SELECT END DO END SUBROUTINE output_run_details SUBROUTINE givens_test(n,ntimes,times,passed) ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision USE givens_rotations, ONLY : d_rot, d_rotm ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0_wp ! .. ! .. Non-Generic Interface Blocks .. INTERFACE FUNCTION elaptime(old_time) ! .. ! .. Function Return Value .. REAL :: elaptime ! .. ! .. Scalar Arguments .. REAL, INTENT (INOUT) :: old_time ! .. END FUNCTION elaptime END INTERFACE INTERFACE SUBROUTINE compute_sol(x,a) ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Array Arguments .. REAL (wp), INTENT (IN) :: a(:,:) REAL (wp), INTENT (INOUT) :: x(:) ! .. END SUBROUTINE compute_sol END INTERFACE ! .. ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: n, ntimes INTEGER, PARAMETER:: mtests=4 ! .. ! .. Array Arguments .. REAL (wp), INTENT (OUT) :: times(mtests) LOGICAL, INTENT (OUT) :: passed(mtests) ! .. ! .. Local Scalars .. REAL (wp) :: c, epsval, s, tmg, tsg REAL :: old_time INTEGER :: i, ierror, itest, j, k, m ! .. ! .. Local Arrays .. REAL (wp), ALLOCATABLE :: a(:,:), b(:,:), d(:), x(:), y(:) REAL (wp) :: param(5) ! .. ! .. External Subroutines .. EXTERNAL drot, drotg, drotm, drotmg ! .. ! .. Intrinsic Functions .. INTRINSIC abs, epsilon, maxval, random_number, sqrt, transpose ! .. epsval = sqrt(epsilon(one)) m = 2*n itest = 1 old_time = 0.0 ! Time using unit strides: ALLOCATE (a(n+1,m),b(n+1,m),d(m),x(n),y(n),STAT=ierror) IF (ierror>0) THEN WRITE (4,'(" Error attempting to allocate arrays for least squares")') STOP END IF CALL random_number(a) ! Construct and apply standard transformations: tsg = elaptime(old_time) DO k = 1, ntimes b = a DO j = 1, n DO i = j + 1, m CALL d_rot(b(j,j),b(j,i),n-j+1,b(j+1,j),1,b(j+1,i),1,c,s) END DO END DO END DO tsg = elaptime(old_time) times(itest) = tsg ! Check solutions: y(1:n) = b(n+1,1:n) CALL compute_sol(y(1:n),transpose(b(1:n,1:n))) tmg = elaptime(old_time) b = a DO k = 1, ntimes a = b d = one DO j = 1, n DO i = j + 1, m ! Construct and apply modified transformation: CALL d_rotm(d(j),d(i),a(j,j),a(j,i),n-j+1,a(j+1,j),1,a(j+1,i),1, & param) END DO END DO END DO tmg = elaptime(old_time) times(itest+1) = tmg ! Check solutions: x(1:n) = a(n+1,1:n) CALL compute_sol(x(1:n),transpose(a(1:n,1:n))) y = (x-y)/maxval(abs(y)) passed(itest) = maxval(abs(y)) <= epsval passed(itest+1) = passed(itest) itest = itest + 2 ! Obtain times with Level-1 BLAS, 2 calls per elimination. ! Construct and apply standard transformations: a = b tsg = elaptime(old_time) DO k = 1, ntimes b = a DO j = 1, n DO i = j + 1, m CALL drotg(b(j,j),b(j,i),c,s) CALL drot(n-j+1,b(j+1,j),1,b(j+1,i),1,c,s) END DO END DO END DO tsg = elaptime(old_time) times(itest) = tsg ! Check solutions: y(1:n) = b(n+1,1:n) CALL compute_sol(y(1:n),transpose(b(1:n,1:n))) tmg = elaptime(old_time) b = a DO k = 1, ntimes a = b d = one DO j = 1, n DO i = j + 1, m ! Construct and apply modified transformation: CALL drotmg(d(j),d(i),a(j,j),a(j,i),param) CALL drotm(n-j+1,a(j+1,j),1,a(j+1,i),1,param) END DO END DO END DO tmg = elaptime(old_time) times(itest+1) = tmg ! Check solutions: x(1:n) = a(n+1,1:n) CALL compute_sol(x(1:n),transpose(a(1:n,1:n))) y = (x-y)/maxval(abs(y)) passed(itest) = maxval(abs(y)) <= epsval passed(itest+1) = passed(itest) itest = itest + 2 DEALLOCATE (a,b,d,x,y, stat=ierror) END SUBROUTINE givens_test SHAR_EOF fi # end of overwriting check if test -f 'Dbenchmark.res' then echo shar: will not over-write existing file "'Dbenchmark.res'" else cat << "SHAR_EOF" > 'Dbenchmark.res' Givens Transformation Benchmark Program Double Precision Version Run started: 18:24 on 10 Oct 2002 Platform id: sun f95 -fast Least squares problem size: 20 by 10 Problem solved 100 times per run Timings averaged over 10 complete runs Timing Details -------------- New MG New SG Orig MG Orig SG Average Times: 0.01 0.01 0.01 0.02 Standard Deviations: 0.00 0.00 0.00 0.00 Normalized: 1.09 1.00 1.36 1.61 Max Times: 0.01 0.01 0.01 0.02 Min Times: 0.01 0.01 0.01 0.02 Normalized Min Times: 1.10 1.00 1.37 1.61 *** Optimal timings obtained using new SG SHAR_EOF fi # end of overwriting check if test -f 'Givens_Rotations.f90' then echo shar: will not over-write existing file "'Givens_Rotations.f90'" else cat << "SHAR_EOF" > 'Givens_Rotations.f90' MODULE givens_rotations ! .. Use Statements .. USE import_kinds, ONLY : sp => single_precision, dp => double_precision ! .. ! .. Generic Interface Blocks .. INTERFACE rot MODULE PROCEDURE s_rot MODULE PROCEDURE d_rot END INTERFACE INTERFACE rotm MODULE PROCEDURE s_rotm MODULE PROCEDURE d_rotm END INTERFACE ! .. CONTAINS ! Modified by R. Hanson 24 April 2002. Use more stable downdating. ! Ref. is from Bjorck, Stewart and Stewart, and Chambers. ! !!!!!!!!!!!!!!!!!!!!!!!!! Documentation !!!!!!!!!!!!!!!!!!!!! ! This module contains a set of routines that efficiently implement ! the calculation and application of a Givens plane rotation using ! both the straightforward and modified methods. ! These routines are intended to replace the blas 1 routines ! drot, drotm, drotg, drotmg, srot, srotm, srotg and srotmg. ! Six user callable routines have been provided. ! Calls to d_rot are intended to be direct replacements for a pair of ! calls to the blas level 1 routines drotg and drot; i.e., the ! calculation of the rotation parameters and the application of the ! rotation to a pair of rows has been amalgamated in a single call. ! s_rot replaces calls to srotg and srot. ! In a similar fashion d_rotm and s_rotm replace pairs of calls to drotg ! and drotg, and srotm and srotmg for the modified version of the ! algorithm. ! The routine *_drot and *_drotm mirror the original blas routines in ! their calling sequences. The data to be transformed is passed to the ! routines as a pair of one dimensional, assumed size arrays and the incx ! and incy stride parameters are also available. These routines may be ! used in Fortran 77 codes that use an array element to define the start ! of an array. ! In addition we have provided a pair of generic routines rot and rotm ! that provide the functionality of the *_rot and *_rotm routines ! but allow the same name to be used for either double or single ! precision data. These routines use assumed size arrays to store ! the row data. ! Experiments were conducted using assumed shape arrays designed to ! shorten the calling sequences by allowing the effect of the use of ! strides to be obtained by using a Fortran 90 array slice; for example ! x(1:1+(k-1)*incx:incx). Also the use of assumed shape arrays would ! allow the length of the data matrices to be determined internally using ! the size function. Unfortunately these codes proved to be extrememly ! inefficient at run time and were abandoned. ! The generic routines are implemented as wrappers to the ! *_rot and *_rotm routines. ! In line with the original blas routines both single and double ! precision versions of the routines are available. We have used a ! consistent naming convention for arguments throughout the routines. ! !!!!!!!!!!!!!!!!!!!!!!!!! Routine Descriptions !!!!!!!!!!!!!!!!!!!!!!! ! S_ROT and D_ROT ! --------------- ! construct a Givens transformation, i.e, they compute ! c = cos(theta) and s = sin(theta) such that either ! ( c s) ( w1 ) ( r ) ! ( ) ( ) = ( ) ! (-s c) ( z1 ) ( 0 ) ! or, when dropping data, they compute the hyperbolic ! transformation matrix of the form ! ( c is) ! (-is c) ! where i is the imaginary unit and this unit is implicitly applied. ! When data is dropped the condition, ! w1**2-z1**2 = (w1-z1)*(w1+z1) > 0 ! is required. ! These routines then apply the computed transformation to the data ! ( x_i ) ! ( ) ! ( y_i ) ! where the ith element of the vector x is ! 1 + (i-1)*incx, if incx >= 0, ! and similarly for y and incy. ! Note that these routines do not allow negative values of incx ! and incy, nor do they compute the z value returned by *rotg ! In addition these routines also allow for row removal which ! was only available with the modified Givens routines *ROTMG ! in the original Blas level 1 library. ! S_ROTM and D_ROTM ! ----------------- ! construct the modified Givens rotation matrix H which ! zeros the second component of the 2-vector ! (w1,z1) = (sqrt(d1)*x1, sqrt(d2)*x2)' ! where d1 and d2 are the reciprocals of the scaling factors ! described in the introduction to the paper and appearing ! in the original blas level 1 routines. ! Row removal mathematically requires d2<0 but this is flagged ! in our code by setting d1<0 for efficiency reasons. ! When data is dropped the condition, ! w1**2-z1**2 = (w1-z1)*(w1+z1) > 0 ! is required. ! and apply it to the 2 x N ! ( x' ) ! matrix ( ), where the elements of the x-vector used are ! ( y' ) ! 1 + (i-1)*incx, if incx >= 0, ! and similarly for y and incy. ! If the value of ! NOTE that these routines require the reciprocal of the d values where ! the original blas routines used the d values directly. Also these ! routines do not allow negative values of incx and incy ! ROT and ROTM ! ------------ ! are Fortran 90 generic interfaces to the routines *_ROT and *_ROTM ! respectively. These routines allow the same routine name to be used ! for either single or double precision data. ! In addition, ROTM uses a structure to return the details of the ! transformation matrix used rather than a floating point array. ! !!!!!!!!!!!!!!!!!!!!!!!!! Parameter Descriptions !!!!!!!!!!!!!!!!!!!!!!! ! ! *_ROT, *_ROTM, ROT and ROTM ! --------------------------- ! (In what follows REAL should be replaced by DOUBLE PRECISION ! when calling the D_ROT and D_ROTM routines. Either REAL or DOUBLE ! PRECISION data may be used with the ROT and ROTM routines provided ! that all the real parameters are of the same type in each call) ! W1,Z1 - REAL (ROT only) ! On entry, these define the two vector whose second element is ! to be annihilated. ! On exit, W1 is set to the value of r (see above) and Z1 is ! set to zero. ! K - INTEGER ! On entry, specifies the number of elements (columns) ! in the two data vectors x and y. ! Unchanged on exit. ! X(*),Y(*) - REAL assumed size arrays ! On entry, these contain the row data to be transformed. ! On exit, unless an error is detected in the input data, ! these contain the transformed data. ! INCX - INTEGER ! On entry, INCX specifies the increment parameter used to step ! through the array SX. Must be positive. ! Unchanged on exit. ! INCY - INTEGER ! On entry, INCY specifies the increment parameter used to step ! through the array SY. Must be positive. ! Unchanged on exit. ! C, S - REAL (ROT only) ! On exit, these are elements of the plane rotation or ! hyperbolic matrix. ! X1, X2 - REAL (ROTM only) ! On entry, contains the values that, when combined with the ! scaling factors RD1 and RD2 (below) give the two vector ! (W1, Z1) whose second element is to be annihilated. ! On exit, X1 contains the transformed value and X2 is ! set to zero. ! RD1, RD2 - REAL (ROTM only) ! On entry, contain the scaling values used for the modified ! transformation. These values are the reciprocal of the values ! (D1 and D2) used in the original level 1 blas routines. ! On exit, contain the scaling values for the transformed vector. ! PARAM(5) - REAL array. (ROTM only) ! On exit, the array param is used to control the form of the ! transformation matrix returned by the modified givens ! algorithm by setting param(1) = sflag where ! sflag = -1.0e0 sflag = 0.0e0 sflag = 1.0e0 ! (sh11 sh12) (1.e0 sh12) ( sh11 1.e0) ! H = ( ) ( ) ( ) ! (sh21 sh22) (sh21 1.e0) (-1.e0 sh22) ! In addition if param(1) = -2 then H is set to I. ! Elements 2--5 of param are used to specify the values of ! sh11, sh21, sh12, sh22 respectively. Values of 1.0e0, -1.0e0, ! or 0.0e0, implied by the value of param(1), are stored ! in param but are not used as multipliers in the transformation ! stage. SUBROUTINE s_rotm(rd1,rd2,x1,y1,k,x,incx,y,incy,param) ! Routine to implement the modified Givens transformation ! This routine is equivalent to the two level 1 blas routines ! srotm and srotmg. ! .. Parameters .. REAL (sp), PARAMETER :: one = 1E0_sp REAL (sp), PARAMETER :: quarter = 25E-2_sp REAL (sp), PARAMETER :: zero = 0E0_sp INTEGER, PARAMETER :: clts = 1, error = -1, rescale = 2, slec = 0 ! .. ! .. Scalar Arguments .. REAL (sp), INTENT (INOUT) :: rd1, rd2, x1, y1 INTEGER, INTENT (IN) :: incx, incy, k ! .. ! .. Array Arguments .. REAL (sp), INTENT (OUT) :: param(5) REAL (sp), INTENT (INOUT) :: x(*), y(*) ! .. ! .. Local Scalars .. REAL (sp) :: gam, gamsq, h11, h12, h21, h22, p1, p2, rgam, sx1, sy1, u INTEGER :: i, kx, ky, n ! .. ! .. Intrinsic Functions .. INTRINSIC abs, huge, min, sqrt, tiny ! .. ! Check input parameters; quick return if illegal IF (rd1==zero .OR. rd2<=zero .OR. (k>0 .AND. (incx<=0 .OR. incy<= & 0))) THEN param(1) = error RETURN END IF n = k kx = 1 ky = 1 sx1 = x1 sy1 = y1 h11 = one h21 = zero h12 = zero h22 = one p2 = rd1*sy1 p1 = rd2*sx1 ! This value is dependent on the underlying ! floating-point arithmetic and need be computed once. ! To make the code thread-safe it is done every time. ! When rescaling is needed gam and rgam are computed. ! gamsq = sqrt(min(huge(one),one/tiny(one))*quarter) gamsq = 4.611686E+18 ! If |c| >= |s|: IF (p1*sx1>=abs(p2*sy1)) THEN param(1) = slec ! There is no transformation necessary. IF (p2==zero) GO TO 10 h21 = -sy1/sx1 h12 = p2/p1 u = one - h12*h21 sx1 = sx1*u rd1 = rd1*u rd2 = rd2*u ! Check if rescaling is required on either diagonal factor. IF (abs(rd1)+rd2>gamsq) GO TO 20 ! Transformation : (1 h12) ! (h21 1) param(2:5) = (/h11, h21, h12, h22/) ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n p1 = x(i) + param(4)*y(i) y(i) = y(i) + param(3)*x(i) x(i) = p1 END DO ELSE ! Non-unit but equal strides: DO i = 1, n p1 = x(kx) + param(4)*y(kx) y(kx) = y(kx) + param(3)*x(kx) x(kx) = p1 kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n p1 = x(kx) + param(4)*y(ky) y(ky) = y(ky) + param(3)*x(kx) x(kx) = p1 kx = kx + incx ky = ky + incy END DO END IF GO TO 10 ELSE h22 = sx1/sy1 h11 = p1/p2 h12 = one h21 = -one u = one + h11*h22 p1 = rd2*u rd2 = rd1*u sx1 = sy1*u rd1 = p1 param(1) = clts ! Check if rescaling is required on either diagonal factor. IF (abs(rd1)+rd2>gamsq) GO TO 20 ! Transformation : ( h11 1) ! (-1 h22) param(2:5) = (/h11, h21, h12, h22/) ! Perhaps only an interchange and sign change is needed. IF (h11==zero .AND. h22==zero) THEN DO i = 1, n p1 = y(ky) y(ky) = -x(kx) x(kx) = p1 kx = kx + incx ky = ky + incy END DO GO TO 10 END IF ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n p1 = y(i) + param(2)*x(i) y(i) = -x(i) + param(5)*y(i) x(i) = p1 END DO ELSE ! Non-unit but equal strides: DO i = 1, n p1 = y(kx) + param(2)*x(kx) y(kx) = -x(kx) + param(5)*y(kx) x(kx) = p1 kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n p1 = y(ky) + param(2)*x(kx) y(ky) = -x(kx) + param(5)*y(ky) x(kx) = p1 kx = kx + incx ky = ky + incy END DO END IF END IF 10 CONTINUE ! The form of the matrix is defined in param(1) when here. param(2:5) = (/h11, h21, h12, h22/) y1 = zero x1 = sx1 RETURN 20 CONTINUE ! Compute gam and rgam if they may be needed for rescaling. gam = sqrt(gamsq) rgam = one/gam IF (abs(rd1)>gamsq) THEN rd1 = (rd1*rgam)*rgam sx1 = sx1*gam h11 = h11*rgam h12 = h12*rgam param(1) = rescale END IF IF (rd2>gamsq) THEN rd2 = (rd2*rgam)*rgam h21 = h21*rgam h22 = h22*rgam param(1) = rescale END IF ! Apply scaled transformation. This is a rare event, after start. DO i = 1, n p1 = h11*x((i-1)*incx+1) + h12*y((i-1)*incy+1) y((i-1)*incy+1) = h21*x((i-1)*incx+1) + h22*y((i-1)*incy+1) x((i-1)*incx+1) = p1 END DO param(2:5) = (/h11, h21, h12, h22/) GO TO 10 END SUBROUTINE s_rotm SUBROUTINE d_rotm(rd1,rd2,x1,y1,k,x,incx,y,incy,param) ! Routine to implement the modified Givens transformation ! This routine is equivalent to the two level 1 blas routines ! drotm and drotmg. ! .. Parameters .. REAL (dp), PARAMETER :: one = 1E0_dp REAL (dp), PARAMETER :: quarter = 25E-2_dp REAL (dp), PARAMETER :: zero = 0E0_dp INTEGER, PARAMETER :: clts = 1, error = -1, rescale = 2, slec = 0 ! .. ! .. Scalar Arguments .. REAL (dp), INTENT (INOUT) :: rd1, rd2, x1, y1 INTEGER, INTENT (IN) :: incx, incy, k ! .. ! .. Array Arguments .. REAL (dp), INTENT (OUT) :: param(5) REAL (dp), INTENT (INOUT) :: x(*), y(*) ! .. ! .. Local Scalars .. REAL (dp) :: gam, gamsq, h11, h12, h21, h22, p1, p2, rgam, sx1, sy1, u INTEGER :: i, kx, ky, n ! .. ! .. Intrinsic Functions .. INTRINSIC abs, huge, min, sqrt, tiny ! .. ! Check input parameters; quick return if illegal IF (rd1==zero .OR. rd2<=zero .OR. (k>0 .AND. (incx<=0 .OR. incy<= & 0))) THEN param(1) = error RETURN END IF n = k kx = 1 ky = 1 sx1 = x1 sy1 = y1 h11 = one h21 = zero h12 = zero h22 = one p2 = rd1*sy1 p1 = rd2*sx1 ! Set the value of gamsq on first call to the ! routine. This value is dependent on the underlying ! floating-point arithmetic and need be computed once. ! To make the code thread-safe it is done every time. ! When rescaling is needed gam and rgam are computed. ! gamsq = sqrt(min(huge(one),one/tiny(one))*quarter) gamsq = 3.3519519824856493D+153 ! If |c| >= |s|: IF (p1*sx1>=abs(p2*sy1)) THEN param(1) = slec ! There is no transformation necessary. IF (p2==zero) GO TO 10 h21 = -sy1/sx1 h12 = p2/p1 u = one - h12*h21 sx1 = sx1*u rd1 = rd1*u rd2 = rd2*u ! Check if rescaling is required on either diagonal factor. IF (abs(rd1)+rd2>gamsq) GO TO 20 ! Transformation : (1 h12) ! (h21 1) param(2:5) = (/h11, h21, h12, h22/) ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n p1 = x(i) + param(4)*y(i) y(i) = y(i) + param(3)*x(i) x(i) = p1 END DO ELSE ! Non-unit but equal strides: DO i = 1, n p1 = x(kx) + param(4)*y(kx) y(kx) = y(kx) + param(3)*x(kx) x(kx) = p1 kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n p1 = x(kx) + param(4)*y(ky) y(ky) = y(ky) + param(3)*x(kx) x(kx) = p1 kx = kx + incx ky = ky + incy END DO END IF GO TO 10 ELSE ! exit with an error return if rd1 < 0 at this point IF (rd1gamsq) GO TO 20 ! Transformation : ( h11 1) ! (-1 h22) param(2:5) = (/h11, h21, h12, h22/) ! Perhaps only an interchange and sign change is needed. IF (h11==zero .AND. h22==zero) THEN DO i = 1, n p1 = y(ky) y(ky) = -x(kx) x(kx) = p1 kx = kx + incx ky = ky + incy END DO GO TO 10 END IF ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n p1 = y(i) + param(2)*x(i) y(i) = -x(i) + param(5)*y(i) x(i) = p1 END DO ELSE ! Non-unit but equal strides: DO i = 1, n p1 = y(kx) + param(2)*x(kx) y(kx) = -x(kx) + param(5)*y(kx) x(kx) = p1 kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n p1 = y(ky) + param(2)*x(kx) y(ky) = -x(kx) + param(5)*y(ky) x(kx) = p1 kx = kx + incx ky = ky + incy END DO END IF END IF 10 CONTINUE ! The form of the matrix is defined in param(1) when here. param(2:5) = (/h11, h21, h12, h22/) y1 = zero x1 = sx1 RETURN 20 CONTINUE ! Compute gam and rgam if they may be needed for rescaling. gam = sqrt(gamsq) rgam = one/gam IF (abs(rd1)>gamsq) THEN rd1 = (rd1*rgam)*rgam sx1 = sx1*gam h11 = h11*rgam h12 = h12*rgam param(1) = rescale END IF IF (rd2>gamsq) THEN rd2 = (rd2*rgam)*rgam h21 = h21*rgam h22 = h22*rgam param(1) = rescale END IF ! Apply scaled transformation. This is a rare event, after start. DO i = 1, n p1 = h11*x((i-1)*incx+1) + h12*y((i-1)*incy+1) y((i-1)*incy+1) = h21*x((i-1)*incx+1) + h22*y((i-1)*incy+1) x((i-1)*incx+1) = p1 END DO param(2:5) = (/h11, h21, h12, h22/) GO TO 10 END SUBROUTINE d_rotm SUBROUTINE s_rot(w1,z1,k,x,incx,y,incy,c,s) ! Routine to implement the standard Givens transformation ! This routine is equivalent to the two level 1 blas routines ! srot and srotg. ! .. Parameters .. REAL (sp), PARAMETER :: half = 5E-1_sp REAL (sp), PARAMETER :: one = 1E0_sp REAL (sp), PARAMETER :: quarter = 25E-2_sp REAL (sp), PARAMETER :: zero = 0E0_sp ! .. ! .. Scalar Arguments .. REAL (sp), INTENT (OUT) :: c, s REAL (sp), INTENT (INOUT) :: w1, z1 INTEGER, INTENT (IN) :: incx, incy, k ! .. ! .. Array Arguments .. REAL (sp), INTENT (INOUT) :: x(*), y(*) ! .. ! .. Local Scalars .. REAL (sp) :: a, b, u, v INTEGER :: i, kx, ky, n ! .. ! .. Intrinsic Functions .. INTRINSIC abs, sqrt ! .. ! Test incx and incy are both positive IF (incx<=0 .OR. incy<=0) THEN c = zero s = zero RETURN END IF n = k a = w1 b = z1 ! If transformation is hyperbolic, code will drop data. IF (n<0) GO TO 10 IF (abs(a)>abs(b)) THEN ! Here ABS(w1) > ABS(z1) u = a + a v = b/u ! Note that u and r have the sign of w1 a = sqrt(quarter+v*v) ! Note that c is positive c = half/a s = v*(c+c) a = a*u ELSE ! Here ABS(w1) <= ABS(z1) IF (b/=zero) THEN u = b + b v = a/u ! Note that u and r have the sign of ! z1 (r is immediately stored in a) a = sqrt(quarter+v*v) ! Note that s is positive s = half/a c = v*(s+s) a = a*u ! Here w1 = z1 = 0. ELSE c = one s = zero GO TO 20 END IF END IF ! Possibly B=0, so no arithmetic is needed. IF (c==one .AND. s==zero) GO TO 20 kx = 1 ky = 1 ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n u = c*x(i) + s*y(i) y(i) = -s*x(i) + c*y(i) x(i) = u END DO ELSE ! Non-unit but equal strides: DO i = 1, n u = c*x(kx) + s*y(kx) y(kx) = -s*x(kx) + c*y(kx) x(kx) = u kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n u = c*x(kx) + s*y(ky) y(ky) = -s*x(kx) + c*y(ky) x(kx) = u kx = kx + incx ky = ky + incy END DO END IF GO TO 20 10 CONTINUE n = -(n+1) ! This is an error condition. IF(ABS(a) <= ABS(b)) THEN c = zero s = zero a = zero RETURN ENDIF ! Will have abs(a) > abs(b) here: s=b/a c=(one-s)*(one+s) ! May have c <= 0, even though mathmatically it should be > 0. ! This is also an error condition. IF(c <= zero) THEN c = zero s = zero a = zero RETURN ENDIF c=sqrt(c) a=c*a ! Will need reciprocal of c to apply transformation: b=one/c kx = 1 ky = 1 ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n u = b*(x(i) - s*y(i)) y(i) = -s*u + c*y(i) x(i) = u END DO ELSE ! Non-unit but equal strides: DO i = 1, n u = b*(x(kx) - s*y(kx)) y(kx) = -s*u + c*y(kx) x(kx) = u kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n u = b*(x(kx) - s*y(ky)) y(ky) = -s*u + c*y(ky) x(kx) = u kx = kx + incx ky = ky + incy END DO END IF 20 CONTINUE z1 = zero w1 = a END SUBROUTINE s_rot SUBROUTINE d_rot(w1,z1,k,x,incx,y,incy,c,s) ! Routine to implement the standard Givens transformation ! This routine is equivalent to the two level 1 blas routines ! drot and drotg. ! .. Parameters .. REAL (dp), PARAMETER :: half = 5E-1_dp REAL (dp), PARAMETER :: one = 1E0_dp REAL (dp), PARAMETER :: quarter = 25E-2_dp REAL (dp), PARAMETER :: zero = 0E0_dp ! .. ! .. Scalar Arguments .. REAL (dp), INTENT (OUT) :: c, s REAL (dp), INTENT (INOUT) :: w1, z1 INTEGER, INTENT (IN) :: incx, incy, k ! .. ! .. Array Arguments .. REAL (dp), INTENT (INOUT) :: x(*), y(*) ! .. ! .. Local Scalars .. REAL (dp) :: a, b, u, v INTEGER :: i, kx, ky, n ! .. ! .. Intrinsic Functions .. INTRINSIC abs, sqrt ! .. ! Test incx and incy are both positive IF (incx<=0 .OR. incy<=0) THEN c = zero s = zero RETURN END IF n = k a = w1 b = z1 ! If transformation is hyperbolic, code will drop data. IF (n<0) GO TO 10 IF (abs(a)>abs(b)) THEN ! Here ABS(w1) > ABS(z1) u = a + a v = b/u ! Note that u and r have the sign of W1 a = sqrt(quarter+v*v) ! Note that c is positive c = half/a s = v*(c+c) a = a*u ELSE ! Here ABS(w1) <= ABS(z1) IF (b/=zero) THEN u = b + b v = a/u ! Note that u and r have the sign of ! z1 (r is immediately stored in a) a = sqrt(quarter+v*v) ! Note that s is positive s = half/a c = v*(s+s) a = a*u ! Here W1 = Z1 = 0. ELSE c = one s = zero GO TO 20 END IF END IF ! Possibly B=0, so no arithmetic is needed. IF (c==one .AND. s==zero) GO TO 20 kx = 1 ky = 1 ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n u = c*x(i) + s*y(i) y(i) = -s*x(i) + c*y(i) x(i) = u END DO ELSE ! Non-unit but equal strides: DO i = 1, n u = c*x(kx) + s*y(kx) y(kx) = -s*x(kx) + c*y(kx) x(kx) = u kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n u = c*x(kx) + s*y(ky) y(ky) = -s*x(kx) + c*y(ky) x(kx) = u kx = kx + incx ky = ky + incy END DO END IF GO TO 20 10 CONTINUE n = -(n+1) IF(ABS(a) <= ABS(b)) THEN c = zero s = zero a = zero RETURN ENDIF ! Will have abs(a) > abs(b) here: s=b/a c=(one-s)*(one+s) ! May have c <= 0, even though mathmatically it should be > 0. ! This is also an error condition. IF(c <= zero) THEN c = zero s = zero a = zero RETURN ENDIF c=sqrt(c) a=c*a ! Will need reciprocal of c to apply transformation: b=one/c kx = 1 ky = 1 ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n u = b*(x(i) - s*y(i)) y(i) = -s*u + c*y(i) x(i) = u END DO ELSE ! Non-unit but equal strides: DO i = 1, n u = b*(x(kx) - s*y(kx)) y(kx) = -s*u + c*y(kx) x(kx) = u kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n u = b*(x(kx) - s*y(ky)) y(ky) = -s*u + c*y(ky) x(kx) = u kx = kx + incx ky = ky + incy END DO END IF 20 CONTINUE z1 = zero w1 = a END SUBROUTINE d_rot END MODULE givens_rotations SHAR_EOF fi # end of overwriting check if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' # # Makefile for benchmark program to triangularize 2n x n random matrices # using Givens transformations. # # This makefile builds the binaries for the single and double precision # versions of the benchmark code. # # To run a particular benchmark run # ./dtimer < data # for the double precision benchmark (use stime for the single precision # version) # # The data file should contain the following lines of data # line 1: the problem size N (this tests for the efficiency of the # transformation using vectors of size 1 .. N-1) # line 2: the number of randomly generated problems to be run, NRAND # line 3: the compiler flags used (printed in the results file) # line 4: the number of times the complete test should be run, NTIMES # (NRAND*NTIMES problems are run in total) # line 5: name of the output file for the results # # For example: # ./dtimer < data_dp # ./stimer < data_sp # # will generate two results files, dbenchmark.res and sbenchmark.res # .SUFFIXES: .f90 # F90 should be set to the name of the Fortran 90 compiler to be used F90 = f95 F90 = /opt/NAGWare/bin/f95 # F90OPTS should be set to compiler flags to be used F90OPTS = -C -g -u -xcheck=stkovf -ansi -XlistE F90OPTS = -g90 -gline -dcfuns -nan -save -C # F90LINK and F90LINKOPTS should be set to the relevant linker and # linker options. For most systems these are the same as for the compiler F90LINK = $(F90) F90LINKOPTS = $(F90OPTS) .f90.o: $(F90) $(F90OPTS) -c $*.f90 .f.o: $(F90) $(F90OPTS) -c $*.f STIMEROBJS=Import_Kinds.o Givens_Rotations.o SBenchmark.o srotmg.o blas.o DTIMEROBJS=Import_Kinds.o Givens_Rotations.o DBenchmark.o drotmg.o blas.o all: dtimer stimer stimer: $(STIMEROBJS) $(F90LINK) -o stimer $(F90LINKOPTS) $(STIMEROBJS) dtimer: $(DTIMEROBJS) $(F90LINK) -o dtimer $(F90LINKOPTS) $(DTIMEROBJS) Import_Kinds.o: ../Src/Import_Kinds.f90 $(F90) $(F90OPTS) -c ../Src/Import_Kinds.f90 Givens_Rotations.o: Import_Kinds.o Givens_Rotations.f90 $(F90) $(F90OPTS) -c Givens_Rotations.f90 DBenchmark.o: DBenchmark.f90 $(F90) $(F90OPTS) -c DBenchmark.f90 SBenchmark.o: SBenchmark.f90 $(F90) $(F90OPTS) -c SBenchmark.f90 blas.o: ../../Remark/Drivers/blas.f $(F90) $(F90OPTS) -c ../../Remark/Drivers/blas.f drotmg.o: ../../Remark/Src/drotmg.f $(F90) $(F90OPTS) -c ../../Remark/Src/drotmg.f srotmg.o: ../../Remark/Src/srotmg.f $(F90) $(F90OPTS) -c ../../Remark/Src/srotmg.f clean: rm -rf ./stimer ./dtimer rm -rf *.o *.g90 *.d *.mod ./a.out core SHAR_EOF fi # end of overwriting check if test -f 'SBenchmark.f90' then echo shar: will not over-write existing file "'SBenchmark.f90'" else cat << "SHAR_EOF" > 'SBenchmark.f90' FUNCTION elaptime(old_time) ! Use this timing function in pairs: ! TS = ELAPTIME(TICKS); { Time something}; TS=ELAPTIME(TICKS) ! The first call starts the clock by recording TICKS. ! (Ignore TS the first call.) ! The second call updates TICKS and computes TS, in Seconds. ! This routine allows for 1-clock roll over. ! This is a Fortran 90 standard routine. It may return 0. if the ! underlying processor has no clock. ! .. ! .. Function Return Value .. REAL :: elaptime ! .. ! .. Scalar Arguments .. REAL, INTENT (INOUT) :: old_time ! .. ! .. Local Scalars .. REAL :: newtime ! .. ! .. Intrinsic Functions .. INTRINSIC cpu_time ! .. CALL cpu_time(newtime) elaptime = newtime-old_time old_time=newtime END FUNCTION elaptime PROGRAM main_timer ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: zero = 0.0_wp LOGICAL, PARAMETER :: details = .FALSE. INTEGER, PARAMETER :: mtests=4 character (3), parameter :: month(12) = (/'Jan', 'Feb',& 'Mar','Apr','May','Jun','Jul','Aug','Sep','Oct', & 'Nov', 'Dec'/) ! .. ! .. Local Scalars .. INTEGER :: i, ierror, m, mon, n, ntimes, repeat_timings LOGICAL :: success CHARACTER (64) :: identification, file_name character (8):: date character (11):: time ! .. ! .. Local Arrays .. REAL (wp) :: maxtim(mtests), mintim(mtests) REAL (wp) :: mtime(mtests), sd(mtests) REAL (wp), ALLOCATABLE :: times(:,:) LOGICAL :: passed(mtests) ! .. ! .. External Functions .. LOGICAL, EXTERNAL :: test_pass ! .. ! .. External Subroutines .. EXTERNAL givens_test, output_run_details ! .. ! .. Intrinsic Functions .. INTRINSIC minloc, minval, real, sum, trim ! .. ! Setting details to true will generate detailed timings of ! all methods for each run. WRITE (*,fmt='(/1x,A)') & ' Benchmark Givens Transformations on 2N by N random matrices:' WRITE (*,fmt='(1x,A)',advance='NO') ' ENTER the number of COLUMNS: >' READ (*,*,err=10,end=7) n WRITE (*,fmt='(1x,A)',advance='NO') & ' ENTER the number of times to test: >' READ (*,*,err=10) ntimes WRITE (*,fmt='(1x,A)',advance='NO') & ' ENTER platform identification comment: >' READ (*,fmt='(A64)',err=10) identification WRITE (*,fmt='(1x,A)',advance='NO') & ' ENTER how many times to repeat the whole test: >' READ (*,*) repeat_timings WRITE (*,fmt='(1x,A)',advance='NO') & ' File name for results: >' READ (*,'(a)')file_name open(unit=4, file=file_name, status='unknown', iostat=ierror) IF (ierror/=0)THEN WRITE (*,'(" Error attempting to open results file: ",a)')& & file_name STOP END IF ! ! Print header information call date_and_time(date,time) read(date(5:6),'(i2)')mon WRITE (4,fmt='(15x,''Givens Transformation Benchmark Program'')') WRITE (4,fmt='(22x,''Single Precision Version''/)') WRITE (4,fmt='(15x,''Run started: '',a2,'':'',a2,'' on '', & &a2,'' '',a3,'' '',a4)') time(1:2), time(3:4), & &date(7:8), month(mon), date(1:4) m = 2*n ALLOCATE (times(mtests,repeat_timings),STAT=ierror) IF (ierror>0) THEN WRITE (4,'(" Error attempting to allocate timing array")') STOP END IF success = .TRUE. DO i = 1, repeat_timings CALL givens_test(n,ntimes,times(1:mtests,i),passed) success = success .AND. test_pass(passed) IF (details) THEN CALL output_run_details(times(1:mtests,i),success) END IF END DO ! Compute average time mtime = zero DO i = 1, repeat_timings mtime = mtime + times(:,i) END DO mtime = mtime/real(repeat_timings,wp) ! and standard deviations sd = zero DO i = 1, repeat_timings sd = sd + (times(:,i)-mtime)**2 END DO DO i = 1, mtests maxtim(i) = maxval(times(i,:)) mintim(i) = minval(times(i,:)) ENDDO ! Output synopsis of runs i = sum(minloc(mtime)) WRITE (4,'(15x,"Platform id: ",a/)') trim(identification) WRITE (4,'(15x,"Least squares problem size: ",i4," by ",i4)') m, n WRITE (4,'(15x,"Problem solved ",i5," times per run")')ntimes WRITE (4,'(15x,"Timings averaged over ",i3," complete runs"/)')& &repeat_timings WRITE (4,fmt='(25x,"Timing Details")') WRITE (4,fmt='(25x,"--------------"/)') WRITE (4,fmt='(28x,"New MG",4x,"New SG",3x,"Orig MG",3x,"Orig SG")') 100 format(a22,":",8f10.2) WRITE (4,100)"Average Times", mtime WRITE (4,100)"Standard Deviations", sd WRITE (4,100)"Normalized", mtime/minval(mtime) write (4,'()') WRITE (4,100)"Max Times", maxtim WRITE (4,100)"Min Times", mintim WRITE (4,100)"Normalized Min Times", mintim/minval(mintim) WRITE (4,fmt='(/5x,A)',advance='NO') & & " *** Optimal timings obtained using " SELECT CASE (i) CASE (1) WRITE (4,'("new MG")') CASE (2) WRITE (4,'("new SG")') CASE (3) WRITE (4,'("original Blas-1 MG")') CASE (4) WRITE (4,'("original Blas-1 SG")') END SELECT WRITE (4,'(/)') DEALLOCATE (times) 7 STOP 10 CONTINUE WRITE (4,'(" Error reading input")') END PROGRAM main_timer SUBROUTINE compute_sol(x,a) ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Array Arguments .. REAL (wp), INTENT (IN) :: a(:,:) REAL (wp), INTENT (INOUT) :: x(:) ! .. ! .. Local Scalars .. INTEGER :: i, j, n ! .. ! .. Intrinsic Functions .. INTRINSIC size ! .. n = size(x) DO j = n, 1, -1 x(j) = x(j)/a(j,j) DO i = 1, j - 1 x(i) = x(i) - a(i,j)*x(j) END DO END DO END SUBROUTINE compute_sol FUNCTION test_pass(passed) RESULT (success) ! .. Function Return Value .. LOGICAL :: success ! .. ! .. Parameters .. INTEGER, PARAMETER :: maxval = 4 ! .. ! .. Array Arguments .. LOGICAL, INTENT (IN) :: passed(maxval) ! .. ! .. Local Scalars .. INTEGER :: i LOGICAL :: ltemp ! .. success = .TRUE. DO i = 1, maxval, 2 ltemp = passed(i) .AND. passed(i+1) success = success .AND. ltemp IF ( .NOT. ltemp) THEN SELECT CASE ((i+1)/2) CASE (1) WRITE (4,*) 'Error in test: new MG then SG' CASE (2) WRITE (4,*) 'Error in test: original Blas-1-1 MG then SG' CASE DEFAULT WRITE (4,*) 'Illegal value of i ', i, ' in test_pass' END SELECT END IF END DO END FUNCTION test_pass SUBROUTINE output_run_details(times,success) ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Parameters .. INTEGER, PARAMETER :: maxval = 4 ! .. ! .. Scalar Arguments .. LOGICAL :: success ! .. ! .. Array Arguments .. REAL (wp) :: times(maxval) ! .. ! .. Local Scalars .. INTEGER :: i ! .. IF (success) THEN WRITE (4,'(" All tests successfully ran")') ELSE WRITE (4,'(" IGNORE TIMINGS: One or more tests failed")') END IF DO i = 1, maxval, 2 SELECT CASE ((i+1)/2) CASE (1) WRITE (4,'(1x,A,1x,2F7.2)') 'Timings: new MG then SG:' & , times(i), times(i+1) CASE (2) WRITE (4,'(1x,A,1x,2F7.2)') & &'Timings: original Blas-1-1 MG then SG:',& ×(i), times(i+1) CASE DEFAULT WRITE (4,*) 'Illegal value of i ', i, ' in test_pass' END SELECT END DO END SUBROUTINE output_run_details SUBROUTINE givens_test(n,ntimes,times,passed) ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision USE givens_rotations, ONLY : s_rot, s_rotm ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0_wp ! .. ! .. Non-Generic Interface Blocks .. INTERFACE FUNCTION elaptime(old_time) ! .. ! .. Function Return Value .. REAL :: elaptime ! .. ! .. Scalar Arguments .. REAL, INTENT (INOUT) :: old_time ! .. END FUNCTION elaptime END INTERFACE INTERFACE SUBROUTINE compute_sol(x,a) ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Array Arguments .. REAL (wp), INTENT (IN) :: a(:,:) REAL (wp), INTENT (INOUT) :: x(:) ! .. END SUBROUTINE compute_sol END INTERFACE ! .. ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: n, ntimes INTEGER, PARAMETER:: mtests=4 ! .. ! .. Array Arguments .. REAL (wp), INTENT (OUT) :: times(mtests) LOGICAL, INTENT (OUT) :: passed(mtests) ! .. ! .. Local Scalars .. REAL (wp) :: c, epsval, s, tmg, tsg REAL :: old_time INTEGER :: i, ierror, itest, j, k, m ! .. ! .. Local Arrays .. REAL (wp), ALLOCATABLE :: a(:,:), b(:,:), d(:), x(:), y(:) REAL (wp) :: param(5) ! .. ! .. External Subroutines .. EXTERNAL srot, srotg, srotm, srotmg ! .. ! .. Intrinsic Functions .. INTRINSIC abs, epsilon, maxval, random_number, sqrt, transpose ! .. epsval = sqrt(epsilon(one)) m = 2*n itest = 1 old_time = 0.0 ! Time using unit strides: ALLOCATE (a(n+1,m),b(n+1,m),d(m),x(n),y(n),STAT=ierror) IF (ierror>0) THEN WRITE (4,'(" Error attempting to allocate arrays for least squares")') STOP END IF CALL random_number(a) ! Construct and apply standard transformations: tsg = elaptime(old_time) DO k = 1, ntimes b = a DO j = 1, n DO i = j + 1, m CALL s_rot(b(j,j),b(j,i),n-j+1,b(j+1,j),1,b(j+1,i),1,c,s) END DO END DO END DO tsg = elaptime(old_time) times(itest) = tsg ! Check solutions: y(1:n) = b(n+1,1:n) CALL compute_sol(y(1:n),transpose(b(1:n,1:n))) tmg = elaptime(old_time) b = a DO k = 1, ntimes a = b d = one DO j = 1, n DO i = j + 1, m ! Construct and apply modified transformation: CALL s_rotm(d(j),d(i),a(j,j),a(j,i),n-j+1,a(j+1,j),1,a(j+1,i),1, & param) END DO END DO END DO tmg = elaptime(old_time) times(itest+1) = tmg ! Check solutions: x(1:n) = a(n+1,1:n) CALL compute_sol(x(1:n),transpose(a(1:n,1:n))) y = (x-y)/maxval(abs(y)) passed(itest) = maxval(abs(y)) <= epsval passed(itest+1) = passed(itest) itest = itest + 2 ! Obtain times with Level-1 BLAS, 2 calls per elimination. ! Construct and apply standard transformations: a = b tsg = elaptime(old_time) DO k = 1, ntimes b = a DO j = 1, n DO i = j + 1, m CALL srotg(b(j,j),b(j,i),c,s) CALL srot(n-j+1,b(j+1,j),1,b(j+1,i),1,c,s) END DO END DO END DO tsg = elaptime(old_time) times(itest) = tsg ! Check solutions: y(1:n) = b(n+1,1:n) CALL compute_sol(y(1:n),transpose(b(1:n,1:n))) tmg = elaptime(old_time) b = a DO k = 1, ntimes a = b d = one DO j = 1, n DO i = j + 1, m ! Construct and apply modified transformation: CALL srotmg(d(j),d(i),a(j,j),a(j,i),param) CALL srotm(n-j+1,a(j+1,j),1,a(j+1,i),1,param) END DO END DO END DO tmg = elaptime(old_time) times(itest+1) = tmg ! Check solutions: x(1:n) = a(n+1,1:n) CALL compute_sol(x(1:n),transpose(a(1:n,1:n))) y = (x-y)/maxval(abs(y)) passed(itest) = maxval(abs(y)) <= epsval passed(itest+1) = passed(itest) itest = itest + 2 DEALLOCATE (a,b,d,x,y, stat=ierror) END SUBROUTINE givens_test SHAR_EOF fi # end of overwriting check if test -f 'Sbenchmark.res' then echo shar: will not over-write existing file "'Sbenchmark.res'" else cat << "SHAR_EOF" > 'Sbenchmark.res' Givens Transformation Benchmark Program Single Precision Version Run started: 18:25 on 10 Oct 2002 Platform id: sun f95 -fast Least squares problem size: 20 by 10 Problem solved 100 times per run Timings averaged over 10 complete runs Timing Details -------------- New MG New SG Orig MG Orig SG Average Times: 0.01 0.01 0.01 0.01 Standard Deviations: 0.00 0.00 0.00 0.00 Normalized: 1.13 1.00 1.35 1.67 Max Times: 0.01 0.01 0.01 0.02 Min Times: 0.01 0.01 0.01 0.01 Normalized Min Times: 1.07 1.00 1.36 1.64 *** Optimal timings obtained using new SG SHAR_EOF fi # end of overwriting check if test -f 'data_dp' then echo shar: will not over-write existing file "'data_dp'" else cat << "SHAR_EOF" > 'data_dp' 10 100 sun f95 -fast 10 dbenchmark.res SHAR_EOF fi # end of overwriting check if test -f 'data_sp' then echo shar: will not over-write existing file "'data_sp'" else cat << "SHAR_EOF" > 'data_sp' 1000 10 sun f95 -fast 1 sbenchmark.res SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'DXrot.f90' then echo shar: will not over-write existing file "'DXrot.f90'" else cat << "SHAR_EOF" > 'DXrot.f90' PROGRAM dp_downdate ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: zero = 0.0_wp ! .. ! .. Local Scalars .. INTEGER :: i,n, ntimes ! .. ! .. Local Arrays .. REAL (wp) :: ulbs(2), maxulbs(2), minulbs(2), avulbs(2) ! .. ! .. External Subroutines .. EXTERNAL downdate_test ! .. ! .. n is the problem size: random matrices of size 2*n x n are used ! .. ntimes defines the number of matrices that are used n = 40 ntimes = 10 ! Print header information WRITE (*,fmt='(/15x,''Downdating Test for Givens Transformation ''/)') WRITE (*,fmt='(12x,i5," random matrices of order ",i5, " are used")') & & ntimes, n WRITE (*,fmt='(15x,''Errors given in multiples of ulb''/)') ! initialize maxulbs = zero minulbs = huge(zero) avulbs = zero DO i=1,ntimes CALL downdate_test(n,ulbs) maxulbs(1) = max(maxulbs(1),ulbs(1)) minulbs(1) = min(minulbs(1),ulbs(1)) maxulbs(2) = max(maxulbs(2),ulbs(2)) minulbs(2) = min(minulbs(2),ulbs(2)) avulbs = avulbs + ulbs END DO WRITE (*,'(26x, ''Average'', 9x, ''Max'',11x,''Min'')') WRITE (*,'(3x, ''NEW SG Downdating'',3e14.4)') & & avulbs(1)/real(ntimes), maxulbs(1),minulbs(1) WRITE (*,'(3x, ''NEW MG Downdating'',3e14.4)') & & avulbs(2)/real(ntimes), maxulbs(2),minulbs(2) STOP 10 CONTINUE WRITE (4,'(" Error reading input")') END PROGRAM dp_downdate SUBROUTINE compute_sol(x,a) ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Array Arguments .. REAL (wp), INTENT (IN) :: a(:,:) REAL (wp), INTENT (INOUT) :: x(:) ! .. ! .. Local Scalars .. INTEGER :: i, j, n ! .. ! .. Intrinsic Functions .. INTRINSIC size ! .. n = size(x) DO j = n, 1, -1 x(j) = x(j)/a(j,j) DO i = 1, j - 1 x(i) = x(i) - a(i,j)*x(j) END DO END DO END SUBROUTINE compute_sol SUBROUTINE downdate_test(n,ulbs) ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision USE givens_rotations, ONLY : d_rot, d_rotm ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0_wp ! .. Non-Generic Interface Blocks .. INTERFACE SUBROUTINE compute_sol(x,a) ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Array Arguments .. REAL (wp), INTENT (IN) :: a(:,:) REAL (wp), INTENT (INOUT) :: x(:) ! .. END SUBROUTINE compute_sol END INTERFACE ! .. ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: n ! .. ! .. Array Arguments .. REAL (wp), INTENT (OUT) :: ulbs(2) ! .. ! .. Local Scalars .. REAL (wp) :: c, d1, s INTEGER :: i, ierror, j ! .. ! .. Local Arrays .. REAL (wp), ALLOCATABLE :: a(:,:), b(:,:), d(:), r(:), x(:), y(:) REAL (wp) :: param(5) ! .. ! .. External Subroutines .. EXTERNAL drot, drotg, drotm, drotmg ! .. ! .. Intrinsic Functions .. INTRINSIC abs, maxval, random_number, transpose ! .. ! Time using unit strides: ALLOCATE (a(n+1,n+1),b(n+1,n+1),r(n+1),d(n+1),x(n),y(n),STAT=ierror) IF (ierror>0) THEN WRITE (4,'(" Error attempting to allocate arrays for least squares")') STOP END IF CALL random_number(a) ! Test add then drop accuracy. Both Modified and Standard approaches are used. ! Only unit strides are tested. CALL random_number(a) b = a ! Reduce problem to triangular form: DO j = 1, n DO i = j + 1, n + 1 CALL d_rot(b(j,j),b(j,i),n-j+1,b(j+1,j),1,b(j+1,i),1,c,s) END DO END DO ! Base solution: y(1:n) = b(n+1,1:n) CALL compute_sol(y(1:n),transpose(b(1:n,1:n))) ! Get a new row. Update and then downdate. CALL random_number(r) b(1:n+1,n+1) = r DO j = 1, n ! ADD DATA CALL d_rot(b(j,j),b(j,n+1),n-j+1,b(j+1,j),1,b(j+1,n+1),1,c,s) END DO b(1:n+1,n+1) = r DO j = 1, n ! DROP DATA CALL d_rot(b(j,j),b(j,n+1),-(n-j+1)-1,b(j+1,j),1,b(j+1,n+1),1,c,s) END DO x(1:n) = b(n+1,1:n) CALL compute_sol(x(1:n),transpose(b(1:n,1:n))) x = (x-y)/maxval(abs(y)) ulbs(1) = maxval(abs(x))/epsilon(one) ! Reduce problem to triangular form: d = one DO j = 1, n DO i = j + 1, n + 1 CALL d_rotm(d(j),d(i),a(j,j),a(j,i),n-j+1,a(j+1,j),1,a(j+1,i),1, & param) END DO END DO ! Get a new row. Update and then downdate. a(1:n+1,n+1) = r d(n+1) = one DO j = 1, n ! ADD DATA CALL d_rotm(d(j),d(n+1),a(j,j),a(j,n+1),n-j+1,a(j+1,j),1,a(j+1,n+1),1, & param) END DO a(1:n+1,n+1) = r d(n+1) = one DO j = 1, n d1 = -d(j) ! DROP DATA CALL d_rotm(d1,d(n+1),a(j,j),a(j,n+1),n-j+1,a(j+1,j),1,a(j+1,n+1),1, & param) d(j) = -d1 END DO x(1:n) = a(n+1,1:n) CALL compute_sol(x(1:n),transpose(a(1:n,1:n))) x = (x-y)/maxval(abs(y)) ulbs(2) = maxval(abs(x))/epsilon(one) DEALLOCATE (a,b,d,r,x,y) END SUBROUTINE downdate_test SHAR_EOF fi # end of overwriting check if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' # # Makefile for # 1. battery tests on the updated *rotmg routines # make remark # 2. the downdating test # make downdate # 3. all the battery tests for both the remark and the new versions # make all # .SUFFIXES: .f90 # F90 should be set to the name of the Fortran 90 compiler to be used F90 = f95 F90 = /opt/NAGWare/bin/f95 # F90OPTS should be set to compiler flags to be used F90OPTS = -C -g -u -xcheck=stkovf -ansi -XlistE F90OPTS = -g90 -gline -dcfuns -nan -save -C # F90LINK and F90LINKOPTS should be set to the relevant linker and # linker options. For most systems these are the same as for the compiler F90LINK = $(F90) F90LINKOPTS = $(F90OPTS) .f90.o: $(F90) $(F90OPTS) -c $*.f90 .f.o: $(F90) $(F90OPTS) -c $*.f D_ROTOBJS= Import_Kinds.o Givens_Rotations.o d_rot_battery.o DROTOBJS= Import_Kinds.o blas.o drot_battery.o D_ROTMOBJS= Import_Kinds.o Givens_Rotations.o d_rotm_battery.o DROTMOBJS= Import_Kinds.o drotmg.o blas.o drotm_battery.o DGROTMOBJS= Import_Kinds.o Givens_Rotations.o rotm_battery.f90 DGROTOBJS= Import_Kinds.o Givens_Rotations.o rot_battery.f90 S_ROTOBJS= Import_Kinds.o Givens_Rotations.o s_rot_battery.o SROTOBJS= Import_Kinds.o blas.o srot_battery.o S_ROTMOBJS= Import_Kinds.o Givens_Rotations.o s_rotm_battery.o SROTMOBJS= Import_Kinds.o srotmg.o blas.o srotm_battery.o SGROTMOBJS= Import_Kinds.o Givens_Rotations.o gsrotm_battery.f90 SGROTOBJS= Import_Kinds.o Givens_Rotations.o gsrot_battery.f90 D_ROTDATA= data1 DROTDATA= ../../Remark/Drivers/data2 D_ROTMDATA= data3 DROTMDATA= ../../Remark/Drivers/data4 DGROTMDATA= data3 DGROTDATA= data1 DDATE_OBJS=Import_Kinds.o Givens_Rotations.o DXrot.o drotmg.o blas.o SDATE_OBJS=Import_Kinds.o Givens_Rotations.o SXrot.o srotmg.o blas.o all: res1 res2 res3 res4 res5 res6 res7 res8 res9 res10 res11 res12\ sdate_res ddate_res remark: res4 res5 res10 res11 downdate: sdate_res ddate_res sdate_res: $(SDATE_OBJS) $(F90LINK) -o sdate $(F90LINKOPTS) $(SDATE_OBJS) ./sdate > sdate_res ddate_res: $(DDATE_OBJS) $(F90LINK) -o ddate $(F90LINKOPTS) $(DDATE_OBJS) ./ddate > ddate_res res1: $(DGROTMOBJS) $(DGROTMDATA) $(F90LINK) $(F90LINKOPTS) $(DGROTMOBJS) ./a.out < $(DGROTMDATA) > res1 res2: $(DGROTOBJS) $(DGROTDATA) $(F90LINK) $(F90LINKOPTS) $(DGROTOBJS) ./a.out < $(DGROTDATA) > res2 res3: $(D_ROTOBJS) $(D_ROTDATA) $(F90LINK) $(F90LINKOPTS) $(D_ROTOBJS) ./a.out < $(D_ROTDATA) > res3 res4: $(DROTOBJS) $(DROTDATA) $(F90LINK) $(F90LINKOPTS) $(DROTOBJS) ./a.out < $(DROTDATA) > res4 res5: $(DROTMOBJS) $(DROTMDATA) $(F90LINK) $(F90LINKOPTS) $(DROTMOBJS) ./a.out < $(DROTMDATA) > res5 res6: $(D_ROTMOBJS) $(D_ROTMDATA) $(F90LINK) $(F90LINKOPTS) $(D_ROTMOBJS) ./a.out < $(D_ROTMDATA) > res6 res7: $(SGROTMOBJS) $(DGROTMDATA) $(F90LINK) $(F90LINKOPTS) $(SGROTMOBJS) ./a.out < $(DGROTMDATA) > res7 res8: $(SGROTOBJS) $(DGROTDATA) $(F90LINK) $(F90LINKOPTS) $(SGROTOBJS) ./a.out < $(DGROTDATA) > res8 res9: $(S_ROTOBJS) $(D_ROTDATA) $(F90LINK) $(F90LINKOPTS) $(S_ROTOBJS) ./a.out < $(D_ROTDATA) > res9 res10: $(SROTOBJS) $(DROTDATA) $(F90LINK) $(F90LINKOPTS) $(SROTOBJS) ./a.out < $(DROTDATA) > res10 res11: $(SROTMOBJS) $(DROTMDATA) $(F90LINK) $(F90LINKOPTS) $(SROTMOBJS) ./a.out < $(DROTMDATA) > res11 res12: $(S_ROTMOBJS) $(D_ROTMDATA) $(F90LINK) $(F90LINKOPTS) $(S_ROTMOBJS) ./a.out < $(D_ROTMDATA) > res12 Import_Kinds.o: ../Src/Import_Kinds.f90 $(F90) $(F90OPTS) -c ../Src/Import_Kinds.f90 Givens_Rotations.o: Import_Kinds.o ../Src/Givens_Rotations.f90 $(F90) $(F90OPTS) -c ../Src/Givens_Rotations.f90 DXrot.o: DXrot.f90 $(F90) $(F90OPTS) -c DXrot.f90 SXrot.o: SXrot.f90 $(F90) $(F90OPTS) -c SXrot.f90 gsrotm_battery.o: gsrotm_battery.f90 $(F90) $(F90OPTS) -c gsrotm_battery.f90 rotm_battery.o: rotm_battery.f90 $(F90) $(F90OPTS) -c rotm_battery.f90 rot_battery.o: rot_battery.f90 $(F90) $(F90OPTS) -c rot_battery.f90 d_rot_battery.o: d_rot_battery.f90 $(F90) $(F90OPTS) -c d_rot_battery.f90 s_rot_battery.o: s_rot_battery.f90 $(F90) $(F90OPTS) -c s_rot_battery.f90 d_rotm_battery.o: d_rotm_battery.f90 $(F90) $(F90OPTS) -c d_rotm_battery.f90 s_rotm_battery.o: s_rotm_battery.f90 $(F90) $(F90OPTS) -c s_rotm_battery.f90 gsrot_battery.o: gsrot_battery.f90 $(F90) $(F90OPTS) -c gsrot_battery.f90 blas.o: ../../Remark/Drivers/blas.f $(F90) $(F90OPTS) -c ../../Remark/Drivers/blas.f drot_battery.o: ../../Remark/Drivers/drot_battery.f90 $(F90) $(F90OPTS) -c ../../Remark/Drivers/drot_battery.f90 drotm_battery.o: ../../Remark/Drivers/drotm_battery.f90 $(F90) $(F90OPTS) -c ../../Remark/Drivers/drotm_battery.f90 srot_battery.o: ../../Remark/Drivers/srot_battery.f90 $(F90) $(F90OPTS) -c ../../Remark/Drivers/srot_battery.f90 srotm_battery.o: ../../Remark/Drivers/srotm_battery.f90 $(F90) $(F90OPTS) -c ../../Remark/Drivers/srotm_battery.f90 drotmg.o: ../../Remark/Src/drotmg.f $(F90) $(F90OPTS) -c ../../Remark/Src/drotmg.f srotmg.o: ../../Remark/Src/srotmg.f $(F90) $(F90OPTS) -c ../../Remark/Src/srotmg.f clean: # remove executables and possible debris rm -rf ./a.out core ./sdate ./ddate # remove various module files -- depends on compiler rm -rf *.o *.g90 *.d *.mod # Remove work files for epcf90 compiler rm -rf work.pc work.pcl SHAR_EOF fi # end of overwriting check if test -f 'SXrot.f90' then echo shar: will not over-write existing file "'SXrot.f90'" else cat << "SHAR_EOF" > 'SXrot.f90' PROGRAM sp_downdate ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: zero = 0.0_wp ! .. ! .. Local Scalars .. INTEGER :: i,n, ntimes ! .. ! .. Local Arrays .. REAL (wp) :: ulbs(2), maxulbs(2), minulbs(2), avulbs(2) ! .. ! .. External Subroutines .. EXTERNAL downdate_test ! .. ! .. n is the problem size: random matrices of size 2*n x n are used ! .. ntimes defines the number of matrices that are used n = 40 ntimes = 10 call random_seed ! Print header information WRITE (*,fmt='(/15x,''Downdating Test for Givens Transformation ''/)') WRITE (*,fmt='(12x,i5," random matrices of order ",i5, " are used")') & & ntimes, n WRITE (*,fmt='(15x,''Errors given in multiples of ulb''/)') ! initialize maxulbs = zero minulbs = huge(zero) avulbs = zero DO i=1,ntimes CALL downdate_test(n,ulbs) maxulbs(1) = max(maxulbs(1),ulbs(1)) minulbs(1) = min(minulbs(1),ulbs(1)) maxulbs(2) = max(maxulbs(2),ulbs(2)) minulbs(2) = min(minulbs(2),ulbs(2)) avulbs = avulbs + ulbs END DO WRITE (*,'(26x, ''Average'', 9x, ''Max'',11x,''Min'')') WRITE (*,'(3x, ''NEW SG Downdating'',3e14.4)') & & avulbs(1)/real(ntimes), maxulbs(1),minulbs(1) WRITE (*,'(3x, ''NEW MG Downdating'',3e14.4)') & & avulbs(2)/real(ntimes), maxulbs(2),minulbs(2) STOP 10 CONTINUE WRITE (4,'(" Error reading input")') END PROGRAM sp_downdate SUBROUTINE compute_sol(x,a) ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Array Arguments .. REAL (wp), INTENT (IN) :: a(:,:) REAL (wp), INTENT (INOUT) :: x(:) ! .. ! .. Local Scalars .. INTEGER :: i, j, n ! .. ! .. Intrinsic Functions .. INTRINSIC size ! .. n = size(x) DO j = n, 1, -1 x(j) = x(j)/a(j,j) DO i = 1, j - 1 x(i) = x(i) - a(i,j)*x(j) END DO END DO END SUBROUTINE compute_sol SUBROUTINE downdate_test(n,ulbs) ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision USE givens_rotations, ONLY : s_rot, s_rotm ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0_wp ! .. Non-Generic Interface Blocks .. INTERFACE SUBROUTINE compute_sol(x,a) ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Array Arguments .. REAL (wp), INTENT (IN) :: a(:,:) REAL (wp), INTENT (INOUT) :: x(:) ! .. END SUBROUTINE compute_sol END INTERFACE ! .. ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: n ! .. ! .. Array Arguments .. REAL (wp), INTENT (OUT) :: ulbs(2) ! .. ! .. Local Scalars .. REAL (wp) :: c, d1, s INTEGER :: i, ierror, j ! .. ! .. Local Arrays .. REAL (wp), ALLOCATABLE :: a(:,:), b(:,:), d(:), r(:), x(:), y(:) REAL (wp) :: param(5) ! .. ! .. External Subroutines .. EXTERNAL drot, drotg, drotm, drotmg ! .. ! .. Intrinsic Functions .. INTRINSIC abs, maxval, random_number, transpose ! .. ! Time using unit strides: ALLOCATE (a(n+1,n+1),b(n+1,n+1),r(n+1),d(n+1),x(n),y(n),STAT=ierror) IF (ierror>0) THEN WRITE (4,'(" Error attempting to allocate arrays for least squares")') STOP END IF CALL random_number(a) ! Test add then drop accuracy. Both Modified and Standard approaches are used. ! Only unit strides are tested. CALL random_number(a) b = a ! Reduce problem to triangular form: DO j = 1, n DO i = j + 1, n + 1 CALL s_rot(b(j,j),b(j,i),n-j+1,b(j+1,j),1,b(j+1,i),1,c,s) END DO END DO ! Base solution: y(1:n) = b(n+1,1:n) CALL compute_sol(y(1:n),transpose(b(1:n,1:n))) ! Get a new row. Update and then downdate. CALL random_number(r) b(1:n+1,n+1) = r DO j = 1, n ! ADD DATA CALL s_rot(b(j,j),b(j,n+1),n-j+1,b(j+1,j),1,b(j+1,n+1),1,c,s) END DO b(1:n+1,n+1) = r DO j = 1, n ! DROP DATA CALL s_rot(b(j,j),b(j,n+1),-(n-j+1)-1,b(j+1,j),1,b(j+1,n+1),1,c,s) END DO x(1:n) = b(n+1,1:n) CALL compute_sol(x(1:n),transpose(b(1:n,1:n))) x = (x-y)/maxval(abs(y)) ulbs(1) = maxval(abs(x))/epsilon(one) ! Reduce problem to triangular form: d = one DO j = 1, n DO i = j + 1, n + 1 CALL s_rotm(d(j),d(i),a(j,j),a(j,i),n-j+1,a(j+1,j),1,a(j+1,i),1, & param) END DO END DO ! Get a new row. Update and then downdate. a(1:n+1,n+1) = r d(n+1) = one DO j = 1, n ! ADD DATA CALL s_rotm(d(j),d(n+1),a(j,j),a(j,n+1),n-j+1,a(j+1,j),1,a(j+1,n+1),1, & param) END DO a(1:n+1,n+1) = r d(n+1) = one DO j = 1, n d1 = -d(j) ! DROP DATA CALL s_rotm(d1,d(n+1),a(j,j),a(j,n+1),n-j+1,a(j+1,j),1,a(j+1,n+1),1, & param) d(j) = -d1 END DO x(1:n) = a(n+1,1:n) CALL compute_sol(x(1:n),transpose(a(1:n,1:n))) x = (x-y)/maxval(abs(y)) ulbs(2) = maxval(abs(x))/epsilon(one) DEALLOCATE (a,b,d,r,x,y) END SUBROUTINE downdate_test SHAR_EOF fi # end of overwriting check if test -f 'd_rot_battery.f90' then echo shar: will not over-write existing file "'d_rot_battery.f90'" else cat << "SHAR_EOF" > 'd_rot_battery.f90' PROGRAM d_rot_test ! Test program for rewritten Givens rotation algorithm. ! This routine combines the two original blas 1 routines ! drot and drotg. ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision USE givens_rotations, ONLY : d_rot ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0E0_wp REAL (wp), PARAMETER :: quarter = 0.25E0_wp REAL (wp), PARAMETER :: three = 3.0E0_wp REAL (wp), PARAMETER :: two = 2.0E0_wp REAL (wp), PARAMETER :: zero = 0.0E0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Local Scalars .. REAL (wp) :: c, s INTEGER :: i, incx, incy, k CHARACTER (80) :: title ! .. ! .. Local Arrays .. REAL (wp) :: d(2), x(2), xv(10), yv(10) ! .. ! .. External Subroutines .. EXTERNAL outres ! .. ! .. Intrinsic Functions .. INTRINSIC sqrt ! .. WRITE (6,'(1x,a,/,a)') & ' Test Results for Fortran 77 compatible version of the', & ' standard Givens transformation (double precision)' ! The tests are driven via a data file 10 READ (5,'(a)',end=20) title WRITE (6,'(//">>> ",a)') title READ (5,*) x(1:2) READ (5,*) d(1:2) ! if k > 0 then there is update vector data READ (5,*) k IF (k>0) THEN READ (5,*) incx, incy xv = zero yv = zero IF (incx>0 .AND. incy>0) THEN READ (5,*) (xv(i),i=1,1+(k-1)*incx,incx) READ (5,*) (yv(i),i=1,1+(k-1)*incy,incy) END IF END IF ! If d(1) is negative in the data file we require row ! removal. This is flagged to d_rot by a negative value ! of k. IF (d(1)0 .AND. incy>0) THEN xv(1:1+(k-1)*incx:incx) = sqrt(d(1))*xv(1:1+(k-1)*incx:incx) yv(1:1+(k-1)*incy:incy) = sqrt(d(2))*yv(1:1+(k-1)*incy:incy) END IF WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) CALL d_rot(x(1),x(2),k,xv,incx,yv,incy,c,s) CALL outres(c,s,k,xv,yv,incx,incy) GO TO 10 20 END PROGRAM d_rot_test SUBROUTINE outres(c,s,k,xv,yv,incx,incy) ! Output routine for results ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: zero = 0.0_wp ! .. ! .. Scalar Arguments .. REAL (wp) :: c, s INTEGER :: incx, incy, k ! .. ! .. Array Arguments .. REAL (wp) :: xv(*), yv(*) ! .. ! .. Local Scalars .. INTEGER :: i ! .. ! .. Intrinsic Functions .. INTRINSIC abs ! .. ! print up results WRITE (6,'(/)') WRITE (6,'(" c and s values:", 2e16.8)') c, s ! error condition is c=s=0 ! print message and exit IF (c==zero .AND. s==zero) THEN WRITE (6,'(/" Error detected in d_rot"/)') ELSE ! Fix k value when row removal is being flagged IF (k<0) THEN k = abs(k) - 1 END IF ! k>0 signifies that data transformation has taken place IF (k>0) THEN WRITE (6,'(/"Output vectors x and y: incx = ",i3,", incy = ",i3)') & incx, incy WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) END IF END IF WRITE (6,'(/"--------------------------------------------------")') END SUBROUTINE outres SHAR_EOF fi # end of overwriting check if test -f 'd_rotm_battery.f90' then echo shar: will not over-write existing file "'d_rotm_battery.f90'" else cat << "SHAR_EOF" > 'd_rotm_battery.f90' PROGRAM d_rotmtest ! Test program for rewritten modified Givens rotation algorithm. ! This routine combines the two original blas 1 routines ! drotm and drotmg. ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision USE givens_rotations, ONLY : d_rotm ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0E0_wp REAL (wp), PARAMETER :: quarter = 0.25E0_wp REAL (wp), PARAMETER :: three = 3.0E0_wp REAL (wp), PARAMETER :: two = 2.0E0_wp REAL (wp), PARAMETER :: zero = 0.0E0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Local Scalars .. REAL (wp) :: gamsq INTEGER :: i, incx, incy, k CHARACTER (80) :: title ! .. ! .. Local Arrays .. REAL (wp) :: d(2), param(5), x(2), xv(10), yv(10) ! .. ! .. External Subroutines .. EXTERNAL outres ! .. ! .. Intrinsic Functions .. INTRINSIC huge, min, tiny ! .. WRITE (6,'(1x,a,/,a)') & ' Test Results for Fortran 77 compatible version of the', & ' modified Givens transformation (double precision)' ! The first tests are driven via a data file ! if k > 0 then there is update vector data 10 READ (5,'(a)',end=20) title WRITE (6,'(//">>> ",a)') title READ (5,*) x(1:2) READ (5,*) d(1:2) ! if k > 0 then there is update vector data READ (5,*) k IF (k>0) THEN READ (5,*) incx, incy xv = zero yv = zero READ (5,*) (xv(i),i=1,1+(k-1)*incx,incx) READ (5,*) (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) ! Generate the reciprocals -- we assume that it is the d ! values that are provided. Anyway in most cases these ! start as one so there is only a risk of confusion ! on exit. IF (d(1)/=zero) THEN d(1) = one/d(1) END IF IF (d(2)/=zero) THEN d(2) = one/d(2) END IF CALL d_rotm(d(1),d(2),x(1),x(2),k,xv,incx,yv,incy,param) CALL outres(param,x,d,k,xv,yv,incx,incy) GO TO 10 ! tests for rescaling and tests for large and small data values 20 CONTINUE gamsq = min(huge(one),one/tiny(one))*quarter ! scaling test d(1) = d(2) < 1/g^2 x(1) = one x(2) = -one d(1) = one/(1.01E0_wp*gamsq) d(2) = d(1) xv(1) = -one xv(2) = two yv(1) = one yv(2) = -two k = 2 d = one/d title = 'Scaling test d(1) = d(2) = 1/(1.01*g^2)' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = d(2) = huge/four; type = anti-diagonal k = -1 x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = huge(one)*quarter d = one/d title = 'd(1) = d(2) = huge/four; Type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = d(2) = two*tiny x(1) = one x(2) = -one d(1) = two*tiny(one) d(2) = two*tiny(one) d = one/d title = 'Scaling test d(1) = d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) =huge/four, d(2) = two*tiny x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = two*tiny(one) d = one/d title = 'Scaling test d(1) = huge/four, d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = two*tiny, d(2) = huge/four x(1) = one x(2) = -one d(1) = two*tiny(one) d(2) = huge(one)*quarter d = one/d title = 'Scaling test d(1) = two*tiny, d(2) = huge/four' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = huge/four, d(2) = one; type = anti-diagonal x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = one d = one/d title = 'd(1) = huge/four, d(2) = one; type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = -huge/four, d(2) = two*tiny x(1) = one x(2) = -one d(1) = -huge(one)*quarter d(2) = two*tiny(one) d = one/d title = 'Scaling test d(1) = -huge/four, d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = -two*tiny, d(2) = huge/four, x(1) = 0; type = anti-diagonal x(1) = zero x(2) = -one d(1) = -two*tiny(one) d(2) = huge(one)*quarter d = one/d title = & 'd(1) = -two*tiny, d(2) = huge/four, x(1) = 0; Type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = -two*tiny, d(2) = huge/four, x(2) = 0; type = diagonal x(1) = one x(2) = zero d(1) = -two*tiny(one) d(2) = huge(one)*quarter d = one/d title = 'd(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL d_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) END PROGRAM d_rotmtest SUBROUTINE outres(param,x,d,k,xv,yv,incx,incy) ! Output routine for the returned values ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Scalar Arguments .. INTEGER :: incx, incy, k ! .. ! .. Array Arguments .. REAL (wp) :: d(2), param(5), x(2), xv(*), yv(*) ! .. ! .. Local Scalars .. INTEGER :: i ! .. ! .. Intrinsic Functions .. INTRINSIC abs, int, sqrt ! .. ! print up results WRITE (6,'(/)') SELECT CASE (int(param(1))) CASE (rescaled) WRITE (6,'(" Type -1, Rescaling has taken place")') WRITE (6,'(" Row 1: ", 2e16.8)') param(2), param(4) WRITE (6,'(" Row 2: ", 2e16.8/)') param(3), param(5) WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') one/d(1:2) CASE (sltc) WRITE (6,'(" Type 0, H is unit diagonal")') WRITE (6,'(/" Off-diagonal elements:", 2e16.8)') param(3), param(4) WRITE (6,'( " x-vector: ", 2e16.8)') x(1:2) WRITE (6,'( " d-vector: ", 2e16.8)') one/d(1:2) CASE (clts) WRITE (6,'(" Type 1, H is unit anti-diagonal")') WRITE (6,'(/" Diagonal elements: ", 2e16.8)') param(2), param(5) WRITE (6,'( " x-vector: ", 2e16.8)') x(1:2) WRITE (6,'( " d-vector: ", 2e16.8)') one/d(1:2) CASE (error) WRITE (6,'(" Error condition detected in Rotmg call")') WRITE (6,'(/"--------------------------------------------------")') RETURN END SELECT ! k>0 if data has been transformed IF (k>0) THEN WRITE (6,'(/"Output vectors x and y: incx = ",i3,", incy = ",i3)') & incx, incy WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) WRITE (6,'(/"Output vectors (scaled to match _rot output)")') xv(1:1+(k-1)*incx:incx) = xv(1:1+(k-1)*incx:incx)/sqrt(abs(d(1))) yv(1:1+(k-1)*incy:incy) = yv(1:1+(k-1)*incy:incy)/sqrt(abs(d(2))) WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(/"--------------------------------------------------")') END SUBROUTINE outres SHAR_EOF fi # end of overwriting check if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << "SHAR_EOF" > 'data1' Straight-forward. Equal unit stride. 2.0 3.0 4.0 9.0 4 1 1 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Straight-forward. Equal non-unit stride 2.0 3.0 4.0 9.0 4 2 2 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Straight-forward. Unequal non-unit stride. 2.0 3.0 4.0 9.0 3 2 3 1.0 2.0 3.0 5.0 6.0 7.0 Straight-forward. 3.0 2.0 9.0 4.0 4 2 2 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 d(1) < zero -- drop data. No application step. 3.0 2.0 -1.0 1.0 0 d(1) < zero -- drop data. Equal unit stride. 3.0 2.0 -1.0 1.0 4 1 1 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 d(1) < zero -- drop data. Equal non-unit stride. 3.0 2.0 -1.0 1.0 4 2 2 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 d(1) < zero -- drop data. Unequal non-unit stride. 3.0 2.0 -1.0 1.0 3 2 3 1.0 2.0 3.0 5.0 6.0 7.0 d(1) < zero (z1>w1) -- c=0, s=0. Error condition. 2.0 3.0 -1.0 1.0 0 Input vector is of the form (x1,0). 2.0 0.0 4.0 9.0 1 1 1 1.0 2.0 Input vector is of the form (0,0). 0.0 0.0 4.0 9.0 0 Input vector is of the form (0,x1). 0.0 3.0 4.0 9.0 0 Error: incx <= 0. 2.0 0.0 4.0 9.0 1 0 1 Error: incy <= 0. 2.0 0.0 4.0 9.0 1 1 0 SHAR_EOF fi # end of overwriting check if test -f 'data3' then echo shar: will not over-write existing file "'data3'" else cat << "SHAR_EOF" > 'data3' Straight-forward. Type = anti-diagonal. Equal unit stride. 2.0 3.0 4.0 9.0 4 1 1 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Straight-forward. Type = anti-diagonal. Equal non-unit stride 2.0 3.0 4.0 9.0 4 2 2 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Straight-forward. Type = anti-diagonal. Unequal non-unit stride. 2.0 3.0 4.0 9.0 3 2 3 1.0 2.0 3.0 5.0 6.0 7.0 Straight-forward. Type = unit diagonal 3.0 2.0 9.0 4.0 4 2 2 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Row removal. H=0. Type = diagonal. Equal unit stride. 3.0 2.0 -1.0 1.0 4 1 1 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Row removal. H=0. Type = diagonal. Equal non-unit stride. 3.0 2.0 -1.0 1.0 4 2 2 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Row removal. H=0. Type = diagonal. Unequal non-unit stride. 3.0 2.0 -1.0 1.0 3 2 3 1.0 2.0 3.0 5.0 6.0 7.0 Attempted row removal but z1>w1. Error. 2.0 3.0 -1.0 1.0 4 1 1 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Input vector is of the form (x1,0). H = I. Type = unit diagonal 2.0 0.0 4.0 9.0 1 1 1 1.0 2.0 Input vector is of the form (0,0). H = I. Type = unit diagonal 0.0 0.0 4.0 9.0 1 1 1 1.0 2.0 Input vector is of the form (0,x1). Type = anti-diagonal. Diag elems = 0 0.0 3.0 4.0 9.0 1 1 1 1.0 2.0 Error: 1/d(2) = 0. 2.0 3.0 4.0 0.0 0 Error: 1/d(1) = 0. 2.0 3.0 0.0 9.0 0 Error: 1/d(2) = 0, x(2) = 0. 2.0 0.0 4.0 0.0 0 Error: 1/d(1) = 0, x(1) = 0. 0.0 3.0 0.0 9.0 0 Error: d(2) < 0. 3.0 2.0 9.0 -4.0 0 SHAR_EOF fi # end of overwriting check if test -f 'ddate_res' then echo shar: will not over-write existing file "'ddate_res'" else cat << "SHAR_EOF" > 'ddate_res' Downdating Test for Givens Transformation 10 random matrices of order 40 are used Errors given in multiples of ulb Average Max Min NEW SG Downdating 0.3457E+03 0.2431E+04 0.1404E+02 NEW MG Downdating 0.2041E+03 0.1181E+04 0.1425E+02 SHAR_EOF fi # end of overwriting check if test -f 'gsrot_battery.f90' then echo shar: will not over-write existing file "'gsrot_battery.f90'" else cat << "SHAR_EOF" > 'gsrot_battery.f90' PROGRAM gsrot_test ! Test program for rewritten Givens rotation algorithm. ! This test call the generic Fortran 90 routine. ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision USE givens_rotations, ONLY : rot ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0E0_wp REAL (wp), PARAMETER :: quarter = 0.25E0_wp REAL (wp), PARAMETER :: three = 3.0E0_wp REAL (wp), PARAMETER :: two = 2.0E0_wp REAL (wp), PARAMETER :: zero = 0.0E0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Local Scalars .. REAL (wp) :: c, s INTEGER :: i, incx, incy, k CHARACTER (80) :: title ! .. ! .. Local Arrays .. REAL (wp) :: d(2), x(2), xv(10), yv(10) ! .. ! .. External Subroutines .. EXTERNAL outres ! .. ! .. Intrinsic Functions .. INTRINSIC sqrt ! .. WRITE (6,'(1x,a,/,a)') & ' Test Results for generic Fortran 90 version of the', & ' standard Givens transformation (single precision)' ! The tests are driven via a data file 10 READ (5,'(a)',end=20) title READ (5,*) x(1:2) READ (5,*) d(1:2) ! if k > 0 then there is update vector data READ (5,*) k IF (k>0) THEN READ (5,*) incx, incy ! incx<1 and incy<1 not sensible -- ignore IF (incx<1 .OR. incy<1) GO TO 10 xv = zero yv = zero READ (5,*) (xv(i),i=1,1+(k-1)*incx,incx) READ (5,*) (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(//">>> ",a)') title ! If d(1) is negative in the data file we require row ! removal. This is flagged to d_rot by a negative value ! of k. IF (d(1) < zero) THEN k = -k - 1 d(1) = -d(1) END IF ! calculate the actual values of w1, z1 from the ! separate d and x values. This allows comparison ! of results between standard and modified methods. x = sqrt(d)*x xv(1:1+(k-1)*incx:incx) = sqrt(d(1))*xv(1:1+(k-1)*incx:incx) yv(1:1+(k-1)*incy:incy) = sqrt(d(2))*yv(1:1+(k-1)*incy:incy) WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) CALL rot(x(1),x(2),k,xv,incx,yv,incy,c,s) CALL outres(c,s,k,xv,yv,incx,incy) GO TO 10 20 END PROGRAM gsrot_test SUBROUTINE outres(c,s,k,xv,yv,incx,incy) ! Output routine for results ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: zero = 0.0_wp ! .. ! .. Scalar Arguments .. REAL (wp) :: c, s INTEGER :: incx, incy, k ! .. ! .. Array Arguments .. REAL (wp) :: xv(*), yv(*) ! .. ! .. Local Scalars .. INTEGER :: i ! .. ! print up results WRITE (6,'(/)') WRITE (6,'(" c and s values:", 2e16.8)') c, s ! error condition is c=s=0 ! print message and exit IF (c==zero .AND. s==zero) THEN WRITE (6,'(/" Error detected in d_rot"/)') ELSE ! k>0 signifies that data transformation has taken place IF (k>0) THEN WRITE (6,'(/"Output vectors x and y: incx = ",i3,", incy = ",i3)') & incx, incy WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) END IF END IF WRITE (6,'(/"--------------------------------------------------")') END SUBROUTINE outres SHAR_EOF fi # end of overwriting check if test -f 'gsrotm_battery.f90' then echo shar: will not over-write existing file "'gsrotm_battery.f90'" else cat << "SHAR_EOF" > 'gsrotm_battery.f90' PROGRAM gsrotmtest ! Test program for rewritten modified Givens rotation algorithm. ! The calls are to the generic Fortran 90 version. ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision USE givens_rotations, ONLY : rotm ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0E0_wp REAL (wp), PARAMETER :: quarter = 0.25E0_wp REAL (wp), PARAMETER :: three = 3.0E0_wp REAL (wp), PARAMETER :: two = 2.0E0_wp REAL (wp), PARAMETER :: zero = 0.0E0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Local Scalars .. REAL (wp) :: gamsq INTEGER :: i, incx, incy, k CHARACTER (80) :: title ! .. ! .. Local Arrays .. REAL (wp) :: d(2), param(5), x(2), xv(10), yv(10) ! .. ! .. External Subroutines .. EXTERNAL outres ! .. ! .. Intrinsic Functions .. INTRINSIC huge, min, tiny ! .. WRITE (6,'(1x,a,/,a)') & ' Test Results for generic Fortran 90 version of the', & ' modified Givens transformation (single precision)' ! main loop 10 READ (5,'(a)',end=20) title WRITE (6,'(//">>> ",a)') title READ (5,*) x(1:2) READ (5,*) d(1:2) ! The first tests are driven via a data file ! if k > 0 then there is update vector data READ (5,*) k IF (k>0) THEN READ (5,*) incx, incy xv = zero yv = zero READ (5,*) (xv(i),i=1,1+(k-1)*incx,incx) READ (5,*) (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) ! Generate the reciprocals -- we assume that it is the d ! values that are provided. Anyway in most cases these ! start as one so there is only a risk of confusion ! on exit. IF (d(1)/=zero) THEN d(1) = one/d(1) END IF IF (d(2)/=zero) THEN d(2) = one/d(2) END IF CALL rotm(d(1),d(2),x(1),x(2),k,xv,incx,yv, incy,param) CALL outres(param,x,d,k,xv,yv,incx,incy) GO TO 10 ! tests for rescaling and tests for large and small data values 20 CONTINUE gamsq = min(huge(one),one/tiny(one))*quarter ! scaling test d(1) = d(2) < 1/g^2 x(1) = one x(2) = -one d(1) = one/(two*gamsq) d(2) = d(1) d = one/d k = 2 xv(1) = -one xv(2) = two yv(1) = one yv(2) = -two title = 'Scaling test d(1) = d(2) = 1/(2*g^2)' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = d(2) = huge/four; type = anti-diagonal x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = huge(one)*quarter d = one/d k = -1 title = 'd(1) = d(2) = huge/four; Type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = d(2) = two*tiny x(1) = one x(2) = -one d(1) = two*tiny(one) d(2) = two*tiny(one) d = one/d k = -1 title = 'Scaling test d(1) = d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) =huge/four, d(2) = two*tiny x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = two*tiny(one) d = one/d k = -1 title = 'Scaling test d(1) = huge/four, d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = two*tiny, d(2) = huge/four x(1) = one x(2) = -one d(1) = two*tiny(one) d(2) = huge(one)*quarter d = one/d k = -1 title = 'Scaling test d(1) = two*tiny, d(2) = huge/four' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = huge/four, d(2) = one; type = anti-diagonal x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = one d = one/d k = -1 title = 'd(1) = huge/four, d(2) = one; type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = -huge/four, d(2) = two*tiny x(1) = one x(2) = -one d(1) = -huge(one)*quarter d(2) = two*tiny(one) d = one/d k = -1 title = 'Scaling test d(1) = -huge/four, d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = -two*tiny, d(2) = huge/four, x(1) = 0; type = anti-diagonal x(1) = zero x(2) = -one d(1) = -two*tiny(one) d(2) = huge(one)*quarter d = one/d k = -1 title = & 'd(1) = -two*tiny, d(2) = huge/four, x(1) = 0; Type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = -two*tiny, d(2) = huge/four, x(2) = 0; type = diagonal x(1) = one x(2) = zero d(1) = -two*tiny(one) d(2) = huge(one)*quarter d = one/d k = -1 title = 'd(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) END PROGRAM gsrotmtest SUBROUTINE outres(param,x,d,k,xv,yv,incx,incy) ! Output routine for the returned values ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Scalar Arguments .. INTEGER :: incx, incy, k ! .. ! .. Array Arguments .. REAL (wp) :: d(2), param(5), x(2), xv(*), yv(*) ! .. ! .. Local Scalars .. INTEGER :: i ! .. ! .. Intrinsic Functions .. INTRINSIC abs, int, sqrt ! .. ! print up results WRITE (6,'(/)') SELECT CASE (int(param(1))) CASE (rescaled) WRITE (6,'(" Type -1, Rescaling has taken place")') WRITE (6,'(" Row 1: ", 2e16.8)') param(2), param(4) WRITE (6,'(" Row 2: ", 2e16.8/)') param(3), param(5) WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') one/d(1:2) CASE (sltc) WRITE (6,'(" Type 0, H is unit diagonal")') WRITE (6,'(/" Off-diagonal elements:", 2e16.8)') param(3), param(4) WRITE (6,'( " x-vector: ", 2e16.8)') x(1:2) WRITE (6,'( " d-vector: ", 2e16.8)') one/d(1:2) CASE (clts) WRITE (6,'(" Type 1, H is unit anti-diagonal")') WRITE (6,'(/" Diagonal elements: ", 2e16.8)') param(2), param(5) WRITE (6,'( " x-vector: ", 2e16.8)') x(1:2) WRITE (6,'( " d-vector: ", 2e16.8)') one/d(1:2) CASE (error) WRITE (6,'(" Error condition detected in Rotmg call")') WRITE (6,'(/"--------------------------------------------------")') RETURN END SELECT ! k>0 if data has been transformed IF (k>0) THEN WRITE (6,'(/"Output vectors x and y: incx = ",i3,", incy = ",i3)') & incx, incy WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) WRITE (6,'(/"Output vectors (scaled to match _rot output)")') xv(1:1+(k-1)*incx:incx) = xv(1:1+(k-1)*incx:incx)/sqrt(abs(d(1))) yv(1:1+(k-1)*incy:incy) = yv(1:1+(k-1)*incy:incy)/sqrt(abs(d(2))) WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(/"--------------------------------------------------")') END SUBROUTINE outres SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << "SHAR_EOF" > 'res1' Test Results for generic Fortran 90 version of the modified Givens transformation (double precision) >>> Straight-forward. Type = anti-diagonal. Equal unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.00000000E+00 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 0.91851852E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 0.13333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Equal non-unit stride Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.00000000E+00 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 0.91851852E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 0.13333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Unequal non-unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.00000000E+00 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 -------------------------------------------------- >>> Straight-forward. Type = unit diagonal Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 0.40000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 0.29629630E+00 x-vector: 0.35925926E+01 0.00000000E+00 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.24814815E+01 0.37777778E+01 0.50740741E+01 0.63703704E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463954E+02 y-vector: 0.79197001E+01 0.85289078E+01 0.91381155E+01 0.97473232E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.00000000E+00 d-vector: -0.18000000E+01 0.18000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 -0.13333333E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 -0.17888544E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 0.71554175E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.00000000E+00 d-vector: -0.18000000E+01 0.18000000E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 -0.13333333E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 -0.17888544E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 0.71554175E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.00000000E+00 d-vector: -0.18000000E+01 0.18000000E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 -------------------------------------------------- >>> Attempted row removal but z1>w1. Error. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Input vector is of the form (x1,0). H = I. Type = unit diagonal Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). H = I. Type = unit diagonal Input Data x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,x1). Type = anti-diagonal. Diag elems = 0 Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.30000000E+01 0.00000000E+00 d-vector: 0.90000000E+01 0.40000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: -0.10000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.60000000E+01 y-vector: -0.20000000E+01 -------------------------------------------------- >>> Error: 1/d(2) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(1) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(2) = 0, x(2) = 0. Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(1) = 0, x(1) = 0. Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(2) < 0. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 -0.40000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Scaling test d(1) = d(2) = 1/(2*g^2) Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22471164+308 0.22471164+308 Type -1, Rescaling has taken place Row 1: 0.17272337E-76 -0.17272337E-76 Row 2: 0.17272337E-76 0.17272337E-76 x-vector: 0.11579209E+78 0.00000000E+00 d-vector: 0.74583407-154 0.74583407-154 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.34544674E-76 0.69089348E-76 y-vector: 0.00000000E+00 0.00000000E+00 Output vectors (scaled to match _rot output) x-vector: -0.29833363-153 0.59666726-153 y-vector: 0.00000000E+00 0.00000000E+00 -------------------------------------------------- >>> d(1) = d(2) = huge/four; Type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22250739-307 0.22250739-307 Type 0, H is unit diagonal Off-diagonal elements: 0.10000000E+01 -0.10000000E+01 x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.22471164+308 0.22471164+308 -------------------------------------------------- >>> Scaling test d(1) = d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22471164+308 0.22471164+308 Type -1, Rescaling has taken place Row 1: 0.17272337E-76 -0.17272337E-76 Row 2: 0.17272337E-76 0.17272337E-76 x-vector: 0.11579209E+78 0.00000000E+00 d-vector: 0.74583407-154 0.74583407-154 -------------------------------------------------- >>> Scaling test d(1) = huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22250739-307 0.22471164+308 Type -1, Rescaling has taken place Row 1: 0.10000000E+01 -0.00000000E+00 Row 2: 0.17272337E-76 0.17272337E-76 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.44942328+308 0.14916681-153 -------------------------------------------------- >>> Scaling test d(1) = two*tiny, d(2) = huge/four Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22471164+308 0.22250739-307 Type -1, Rescaling has taken place Row 1: -0.00000000E+00 0.10000000E+01 Row 2: -0.17272337E-76 -0.17272337E-76 x-vector: -0.10000000E+01 0.00000000E+00 d-vector: 0.44942328+308 0.14916681-153 -------------------------------------------------- >>> d(1) = huge/four, d(2) = one; type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22250739-307 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.10000000E+01 -0.22250739-307 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.44942328+308 0.10000000E+01 -------------------------------------------------- >>> Scaling test d(1) = -huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: -0.22250739-307 0.22471164+308 Type -1, Rescaling has taken place Row 1: 0.10000000E+01 0.00000000E+00 Row 2: 0.17272337E-76 0.17272337E-76 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.44942328+308 0.14916681-153 -------------------------------------------------- >>> d(1) = -two*tiny, d(2) = huge/four, x(1) = 0; Type = anti-diagonal Input Data x-vector: 0.00000000E+00 -0.10000000E+01 d-vector: -0.22471164+308 0.22250739-307 Error condition detected in Rotmg call -------------------------------------------------- >>> d(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal Input Data x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.22471164+308 0.22250739-307 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.44501477-307 0.44942328+308 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res10' then echo shar: will not over-write existing file "'res10'" else cat << "SHAR_EOF" > 'res10' Test Results for BLAS version of the standard Givens transformation (single precision) >>> Straight-forward. Equal unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381156E+00 x-vector: 0.98488579E+01 0.24622145E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180584E+02 y-vector: 0.42644539E+01 0.36552458E+01 0.30460386E+01 0.24368305E+01 -------------------------------------------------- >>> Straight-forward. Equal non-unit stride Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381156E+00 x-vector: 0.98488579E+01 0.24622145E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180584E+02 y-vector: 0.42644539E+01 0.36552458E+01 0.30460386E+01 0.24368305E+01 -------------------------------------------------- >>> Straight-forward. Unequal non-unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381156E+00 x-vector: 0.98488579E+01 0.24622145E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552458E+01 0.30460386E+01 -------------------------------------------------- >>> Straight-forward. Input Data x-vector: 0.90000000E+01 0.40000000E+01 c and s values: 0.91381156E+00 0.40613845E+00 x-vector: 0.98488579E+01 0.40613845E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463953E+02 y-vector: 0.79197006E+01 0.85289078E+01 0.91381159E+01 0.97473240E+01 -------------------------------------------------- >>> Input vector is of the form (x1,0). Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 x-vector: 0.40000000E+01 0.00000000E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). Input Data x-vector: 0.00000000E+00 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 x-vector: 0.00000000E+00 0.00000000E+00 -------------------------------------------------- >>> Input vector is of the form (0,x1). Input Data x-vector: 0.00000000E+00 0.90000000E+01 c and s values: 0.00000000E+00 0.10000000E+01 x-vector: 0.90000000E+01 0.10000000E+01 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res11' then echo shar: will not over-write existing file "'res11'" else cat << "SHAR_EOF" > 'res11' Test Results for BLAS version of the modified Givens transformation (single precision) >>> Straight-forward. Type = anti-diagonal. Equal unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.30000000E+01 d-vector: 0.75154638E+01 0.33402061E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 0.91851854E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 0.13333335E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519449E+02 0.18073162E+02 0.21626873E+02 0.25180586E+02 y-vector: 0.42644544E+01 0.36552463E+01 0.30460391E+01 0.24368312E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Equal non-unit stride Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.30000000E+01 d-vector: 0.75154638E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 0.91851854E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 0.13333335E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519449E+02 0.18073162E+02 0.21626873E+02 0.25180586E+02 y-vector: 0.42644544E+01 0.36552463E+01 0.30460391E+01 0.24368312E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Unequal non-unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.30000000E+01 d-vector: 0.75154638E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519449E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644544E+01 0.36552463E+01 0.30460391E+01 -------------------------------------------------- >>> Straight-forward. Type = unit diagonal Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 0.40000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 0.29629630E+00 x-vector: 0.35925927E+01 0.20000000E+01 d-vector: 0.75154638E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.24814816E+01 0.37777777E+01 0.50740738E+01 0.63703704E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: 0.68028193E+01 0.10356530E+02 0.13910241E+02 0.17463953E+02 y-vector: 0.79197006E+01 0.85289078E+01 0.91381159E+01 0.97473230E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.10000000E+01 -0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.20000000E+01 d-vector: 0.18000001E+01 -0.18000001E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 -0.13333335E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 -0.17888546E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.10000000E+01 -0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.20000000E+01 d-vector: 0.18000001E+01 -0.18000001E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 -0.13333335E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 -0.17888546E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.10000000E+01 -0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.20000000E+01 d-vector: 0.18000001E+01 -0.18000001E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 -------------------------------------------------- >>> Attempted row removal but z1>w1. Error. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.10000000E+01 -0.10000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Input vector is of the form (x1,0). H = I. Type = unit diagonal Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type -2, H is unit matrix x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). H = I. Type = unit diagonal Input Data x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type -2, H is unit matrix x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,x1). Type = anti-diagonal. Diag elems = 0 Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.30000000E+01 0.30000000E+01 d-vector: 0.90000000E+01 0.40000000E+01 -------------------------------------------------- >>> Error: d(2) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(1) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(2) = 0. Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(1) = 0. Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Scaling test d(1) = d(2) = 1/(2*g^2) Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.23509887E-37 0.23509887E-37 Type -1, Rescaling has taken place Row 1: -0.21684043E-18 0.21684043E-18 Row 2: -0.21684043E-18 -0.21684043E-18 x-vector: -0.43368087E-18 -0.10000000E+01 d-vector: 0.25000000E+00 0.25000000E+00 -------------------------------------------------- >>> d(1) = d(2) = huge/four; Type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.85070587E+38 0.85070587E+38 Type -1, Rescaling has taken place Row 1: -0.46116860E+19 0.46116860E+19 Row 2: -0.46116860E+19 -0.46116860E+19 x-vector: -0.92233720E+19 -0.10000000E+01 d-vector: 0.19999999E+01 0.19999999E+01 -------------------------------------------------- >>> Scaling test d(1) = d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.42535296E+38 0.42535296E+38 Type -1, Rescaling has taken place Row 1: -0.46116860E+19 0.46116860E+19 Row 2: -0.46116860E+19 -0.46116860E+19 x-vector: -0.92233720E+19 -0.10000000E+01 d-vector: 0.10000000E+01 0.10000000E+01 -------------------------------------------------- >>> Scaling test d(1) = huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.11754945E-37 0.42535296E+38 Type -1, Rescaling has taken place Row 1: -0.00000000E+00 0.46116860E+19 Row 2: -0.21684043E-18 -0.21684043E-18 x-vector: -0.46116860E+19 -0.10000000E+01 d-vector: 0.20000000E+01 0.25000003E+00 -------------------------------------------------- >>> Scaling test d(1) = two*tiny, d(2) = huge/four Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.42535296E+38 0.11754945E-37 Type -1, Rescaling has taken place Row 1: 0.46116860E+19 -0.00000000E+00 Row 2: 0.21684043E-18 0.21684043E-18 x-vector: 0.46116860E+19 -0.10000000E+01 d-vector: 0.20000000E+01 0.25000003E+00 -------------------------------------------------- >>> d(1) = huge/four, d(2) = one; type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.85070587E+38 0.10000000E+01 Type -1, Rescaling has taken place Row 1: 0.46116860E+19 -0.54210115E-19 Row 2: 0.10000000E+01 0.10000000E+01 x-vector: 0.46116860E+19 -0.10000000E+01 d-vector: 0.39999998E+01 0.10000000E+01 -------------------------------------------------- >>> Scaling test d(1) = -huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.85070587E+38 -0.23509887E-37 Type -1, Rescaling has taken place Row 1: 0.46116860E+19 0.00000000E+00 Row 2: 0.21684043E-18 0.21684043E-18 x-vector: 0.46116860E+19 -0.10000000E+01 d-vector: 0.39999998E+01 -0.50000000E+00 -------------------------------------------------- >>> d(1) = two*tiny, d(2) = -huge/four, x(1) = 0; Type = anti-diagonal Input Data x-vector: 0.00000000E+00 -0.10000000E+01 d-vector: 0.23509887E-37 -0.85070587E+38 Error condition detected in Rotmg call -------------------------------------------------- >>> d(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal Input Data x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.23509887E-37 -0.85070587E+38 Type -2, H is unit matrix x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.23509887E-37 -0.85070587E+38 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res12' then echo shar: will not over-write existing file "'res12'" else cat << "SHAR_EOF" > 'res12' Test Results for Fortran 77 compatible version of the modified Givens transformation (single precision) >>> Straight-forward. Type = anti-diagonal. Equal unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.00000000E+00 d-vector: 0.75154643E+01 0.33402061E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 0.91851854E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 0.13333335E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180586E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460391E+01 0.24368310E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Equal non-unit stride Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.00000000E+00 d-vector: 0.75154643E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 0.91851854E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 0.13333335E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180586E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460391E+01 0.24368310E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Unequal non-unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.00000000E+00 d-vector: 0.75154643E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460391E+01 -------------------------------------------------- >>> Straight-forward. Type = unit diagonal Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 0.40000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 0.29629630E+00 x-vector: 0.35925927E+01 0.00000000E+00 d-vector: 0.75154643E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.24814816E+01 0.37777777E+01 0.50740738E+01 0.63703704E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: 0.68028197E+01 0.10356531E+02 0.13910242E+02 0.17463955E+02 y-vector: 0.79197001E+01 0.85289078E+01 0.91381149E+01 0.97473221E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.00000000E+00 d-vector: -0.18000001E+01 0.18000001E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 -0.13333335E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 -0.17888547E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.00000000E+00 d-vector: -0.18000001E+01 0.18000001E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 -0.13333335E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 -0.17888547E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.00000000E+00 d-vector: -0.18000001E+01 0.18000001E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 -------------------------------------------------- >>> Attempted row removal but z1>w1. Error. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: -0.66666669E+00 0.66666669E+00 x-vector: 0.16666665E+01 0.00000000E+00 d-vector: 0.18000001E+01 -0.18000001E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 0.13333335E+01 Output vectors (scaled to match _rot output) x-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 y-vector: 0.31304955E+01 0.26832817E+01 0.22360685E+01 0.17888547E+01 -------------------------------------------------- >>> Input vector is of the form (x1,0). H = I. Type = unit diagonal Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). H = I. Type = unit diagonal Input Data x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,x1). Type = anti-diagonal. Diag elems = 0 Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.30000000E+01 0.00000000E+00 d-vector: 0.90000000E+01 0.40000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: -0.10000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.60000000E+01 y-vector: -0.20000000E+01 -------------------------------------------------- >>> Error: 1/d(2) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(1) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(2) = 0, x(2) = 0. Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(1) = 0, x(1) = 0. Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(2) < 0. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 -0.40000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Scaling test d(1) = d(2) = 1/(1.01E0_wp*g^2) Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.42535296E+38 0.42535296E+38 Type -1, Rescaling has taken place Row 1: 0.46566129E-09 -0.46566129E-09 Row 2: 0.46566129E-09 0.46566129E-09 x-vector: 0.42949673E+10 0.00000000E+00 d-vector: 0.54210109E-19 0.54210109E-19 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.93132257E-09 0.18626451E-08 y-vector: 0.00000000E+00 0.00000000E+00 Output vectors (scaled to match _rot output) x-vector: -0.21684043E-18 0.43368087E-18 y-vector: 0.00000000E+00 0.00000000E+00 -------------------------------------------------- >>> d(1) = d(2) = huge/four; Type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.11754945E-37 0.11754945E-37 Type 0, H is unit diagonal Off-diagonal elements: 0.10000000E+01 -0.10000000E+01 x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.42535291E+38 0.42535291E+38 -------------------------------------------------- >>> Scaling test d(1) = d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.42535296E+38 0.42535296E+38 Type -1, Rescaling has taken place Row 1: 0.46566129E-09 -0.46566129E-09 Row 2: 0.46566129E-09 0.46566129E-09 x-vector: 0.42949673E+10 0.00000000E+00 d-vector: 0.54210109E-19 0.54210109E-19 -------------------------------------------------- >>> Scaling test d(1) = huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.11754945E-37 0.42535296E+38 Type -1, Rescaling has taken place Row 1: 0.10000000E+01 -0.00000000E+00 Row 2: 0.46566129E-09 0.46566129E-09 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.85070582E+38 0.10842022E-18 -------------------------------------------------- >>> Scaling test d(1) = two*tiny, d(2) = huge/four Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.42535296E+38 0.11754945E-37 Type -1, Rescaling has taken place Row 1: -0.00000000E+00 0.10000000E+01 Row 2: -0.46566129E-09 -0.46566129E-09 x-vector: -0.10000000E+01 0.00000000E+00 d-vector: 0.85070582E+38 0.10842022E-18 -------------------------------------------------- >>> d(1) = huge/four, d(2) = one; type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.11754945E-37 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.10000000E+01 -0.11754945E-37 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.85070582E+38 0.10000000E+01 -------------------------------------------------- >>> Scaling test d(1) = -huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: -0.11754945E-37 0.42535296E+38 Type -1, Rescaling has taken place Row 1: 0.10000000E+01 0.00000000E+00 Row 2: 0.46566129E-09 0.46566129E-09 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.85070582E+38 0.10842022E-18 -------------------------------------------------- >>> d(1) = -two*tiny, d(2) = huge/four, x(1) = 0; Type = anti-diagonal Input Data x-vector: 0.00000000E+00 -0.10000000E+01 d-vector: -0.42535296E+38 0.11754945E-37 Type 1, H is unit anti-diagonal Diagonal elements: 0.00000000E+00 -0.00000000E+00 x-vector: -0.10000000E+01 0.00000000E+00 d-vector: 0.85070582E+38 -0.23509887E-37 -------------------------------------------------- >>> d(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal Input Data x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.42535296E+38 0.11754945E-37 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.23509887E-37 0.85070582E+38 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res2' then echo shar: will not over-write existing file "'res2'" else cat << "SHAR_EOF" > 'res2' Test Results for generic Fortran 90 version of the standard Givens transformation (double precision) >>> Straight-forward. Equal unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Equal non-unit stride Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Unequal non-unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 -------------------------------------------------- >>> Straight-forward. Input Data x-vector: 0.90000000E+01 0.40000000E+01 c and s values: 0.91381155E+00 0.40613847E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463954E+02 y-vector: 0.79197001E+01 0.85289078E+01 0.91381155E+01 0.97473232E+01 -------------------------------------------------- >>> d(1) < zero -- drop data. No application step. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535599E+00 0.66666667E+00 -------------------------------------------------- >>> d(1) < zero -- drop data. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535599E+00 0.66666667E+00 -------------------------------------------------- >>> d(1) < zero -- drop data. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535599E+00 0.66666667E+00 -------------------------------------------------- >>> d(1) < zero -- drop data. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535599E+00 0.66666667E+00 -------------------------------------------------- >>> d(1) < zero (z1>w1) -- c=0, s=0. Error condition. Input Data x-vector: 0.20000000E+01 0.30000000E+01 c and s values: 0.00000000E+00 0.00000000E+00 Error detected in d_rot -------------------------------------------------- >>> Input vector is of the form (x1,0). Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). Input Data x-vector: 0.00000000E+00 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 -------------------------------------------------- >>> Input vector is of the form (0,x1). Input Data x-vector: 0.00000000E+00 0.90000000E+01 c and s values: 0.00000000E+00 0.10000000E+01 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res3' then echo shar: will not over-write existing file "'res3'" else cat << "SHAR_EOF" > 'res3' Test Results for Fortran 77 compatible version of the standard Givens transformation (double precision) >>> Straight-forward. Equal unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Equal non-unit stride Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Unequal non-unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 -------------------------------------------------- >>> Straight-forward. Input Data x-vector: 0.90000000E+01 0.40000000E+01 c and s values: 0.91381155E+00 0.40613847E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463954E+02 y-vector: 0.79197001E+01 0.85289078E+01 0.91381155E+01 0.97473232E+01 -------------------------------------------------- >>> d(1) < zero -- drop data. No application step. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535599E+00 0.66666667E+00 -------------------------------------------------- >>> d(1) < zero -- drop data. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535599E+00 0.66666667E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 -0.17888544E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 0.71554175E+01 -------------------------------------------------- >>> d(1) < zero -- drop data. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535599E+00 0.66666667E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 -0.17888544E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 0.71554175E+01 -------------------------------------------------- >>> d(1) < zero -- drop data. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535599E+00 0.66666667E+00 Output vectors x and y: incx = 2, incy = 3 x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 -------------------------------------------------- >>> d(1) < zero (z1>w1) -- c=0, s=0. Error condition. Input Data x-vector: 0.20000000E+01 0.30000000E+01 c and s values: 0.00000000E+00 0.00000000E+00 Error detected in d_rot -------------------------------------------------- >>> Input vector is of the form (x1,0). Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). Input Data x-vector: 0.00000000E+00 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 -------------------------------------------------- >>> Input vector is of the form (0,x1). Input Data x-vector: 0.00000000E+00 0.90000000E+01 c and s values: 0.00000000E+00 0.10000000E+01 -------------------------------------------------- >>> Error: incx <= 0. Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.00000000E+00 0.00000000E+00 Error detected in d_rot -------------------------------------------------- >>> Error: incy <= 0. Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.00000000E+00 0.00000000E+00 Error detected in d_rot -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res4' then echo shar: will not over-write existing file "'res4'" else cat << "SHAR_EOF" > 'res4' Test Results for BLAS version of the standard Givens transformation (double precision) >>> Straight-forward. Equal unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 x-vector: 0.98488578E+01 0.24622145E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Equal non-unit stride Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 x-vector: 0.98488578E+01 0.24622145E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Unequal non-unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613847E+00 0.91381155E+00 x-vector: 0.98488578E+01 0.24622145E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 -------------------------------------------------- >>> Straight-forward. Input Data x-vector: 0.90000000E+01 0.40000000E+01 c and s values: 0.91381155E+00 0.40613847E+00 x-vector: 0.98488578E+01 0.40613847E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463954E+02 y-vector: 0.79197001E+01 0.85289078E+01 0.91381155E+01 0.97473232E+01 -------------------------------------------------- >>> Input vector is of the form (x1,0). Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 x-vector: 0.40000000E+01 0.00000000E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). Input Data x-vector: 0.00000000E+00 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 x-vector: 0.00000000E+00 0.00000000E+00 -------------------------------------------------- >>> Input vector is of the form (0,x1). Input Data x-vector: 0.00000000E+00 0.90000000E+01 c and s values: 0.00000000E+00 0.10000000E+01 x-vector: 0.90000000E+01 0.10000000E+01 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res5' then echo shar: will not over-write existing file "'res5'" else cat << "SHAR_EOF" > 'res5' Test Results for BLAS version of the modified Givens transformation (double precision) >>> Straight-forward. Type = anti-diagonal. Equal unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.30000000E+01 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 0.91851852E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 0.13333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Equal non-unit stride Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.30000000E+01 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 0.91851852E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 0.13333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Unequal non-unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.30000000E+01 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 -------------------------------------------------- >>> Straight-forward. Type = unit diagonal Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 0.40000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 0.29629630E+00 x-vector: 0.35925926E+01 0.20000000E+01 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.24814815E+01 0.37777778E+01 0.50740741E+01 0.63703704E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463954E+02 y-vector: 0.79197001E+01 0.85289078E+01 0.91381155E+01 0.97473232E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.10000000E+01 -0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.20000000E+01 d-vector: 0.18000000E+01 -0.18000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 -0.13333333E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 -0.17888544E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 0.71554175E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.10000000E+01 -0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.20000000E+01 d-vector: 0.18000000E+01 -0.18000000E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 -0.13333333E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 -0.17888544E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 0.71554175E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.10000000E+01 -0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.20000000E+01 d-vector: 0.18000000E+01 -0.18000000E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 -------------------------------------------------- >>> Attempted row removal but z1>w1. Error. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.10000000E+01 -0.10000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Input vector is of the form (x1,0). H = I. Type = unit diagonal Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type -2, H is unit matrix x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). H = I. Type = unit diagonal Input Data x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type -2, H is unit matrix x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,x1). Type = anti-diagonal. Diag elems = 0 Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.30000000E+01 0.30000000E+01 d-vector: 0.90000000E+01 0.40000000E+01 -------------------------------------------------- >>> Error: d(2) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(1) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(2) = 0. Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(1) = 0. Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Scaling test d(1) = d(2) = 1/(2*g^2) Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.44501477-307 0.44501477-307 Type -1, Rescaling has taken place Row 1: -0.29833363-153 0.29833363-153 Row 2: -0.29833363-153 -0.29833363-153 x-vector: -0.59666726-153 -0.10000000E+01 d-vector: 0.25000000E+00 0.25000000E+00 -------------------------------------------------- >>> d(1) = d(2) = huge/four; Type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.44942328+308 0.44942328+308 Type -1, Rescaling has taken place Row 1: -0.33519520+154 0.33519520+154 Row 2: -0.33519520+154 -0.33519520+154 x-vector: -0.67039040+154 -0.10000000E+01 d-vector: 0.20000000E+01 0.20000000E+01 -------------------------------------------------- >>> Scaling test d(1) = d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22471164+308 0.22471164+308 Type -1, Rescaling has taken place Row 1: -0.33519520+154 0.33519520+154 Row 2: -0.33519520+154 -0.33519520+154 x-vector: -0.67039040+154 -0.10000000E+01 d-vector: 0.10000000E+01 0.10000000E+01 -------------------------------------------------- >>> Scaling test d(1) = huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22250739-307 0.22471164+308 Type -1, Rescaling has taken place Row 1: -0.00000000E+00 0.33519520+154 Row 2: -0.29833363-153 -0.29833363-153 x-vector: -0.33519520+154 -0.10000000E+01 d-vector: 0.20000000E+01 0.25000000E+00 -------------------------------------------------- >>> Scaling test d(1) = two*tiny, d(2) = huge/four Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22471164+308 0.22250739-307 Type -1, Rescaling has taken place Row 1: 0.33519520+154 -0.00000000E+00 Row 2: 0.29833363-153 0.29833363-153 x-vector: 0.33519520+154 -0.10000000E+01 d-vector: 0.20000000E+01 0.25000000E+00 -------------------------------------------------- >>> d(1) = huge/four, d(2) = one; type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.44942328+308 0.10000000E+01 Type -1, Rescaling has taken place Row 1: 0.33519520+154 -0.74583407-154 Row 2: 0.10000000E+01 0.10000000E+01 x-vector: 0.33519520+154 -0.10000000E+01 d-vector: 0.40000000E+01 0.10000000E+01 -------------------------------------------------- >>> Scaling test d(1) = -huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.44942328+308 -0.44501477-307 Type -1, Rescaling has taken place Row 1: 0.33519520+154 0.00000000E+00 Row 2: 0.29833363-153 0.29833363-153 x-vector: 0.33519520+154 -0.10000000E+01 d-vector: 0.40000000E+01 -0.50000000E+00 -------------------------------------------------- >>> d(1) = two*tiny, d(2) = -huge/four, x(1) = 0; Type = anti-diagonal Input Data x-vector: 0.00000000E+00 -0.10000000E+01 d-vector: 0.44501477-307 -0.44942328+308 Error condition detected in Rotmg call -------------------------------------------------- >>> d(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal Input Data x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.44501477-307 -0.44942328+308 Type -2, H is unit matrix x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.44501477-307 -0.44942328+308 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res6' then echo shar: will not over-write existing file "'res6'" else cat << "SHAR_EOF" > 'res6' Test Results for Fortran 77 compatible version of the modified Givens transformation (double precision) >>> Straight-forward. Type = anti-diagonal. Equal unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.00000000E+00 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 0.91851852E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 0.13333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Equal non-unit stride Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.00000000E+00 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 0.91851852E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 0.13333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180585E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 0.24368308E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Unequal non-unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666667E+00 x-vector: 0.35925926E+01 0.00000000E+00 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.52962963E+01 0.65925926E+01 0.78888889E+01 y-vector: 0.23333333E+01 0.20000000E+01 0.16666667E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552462E+01 0.30460385E+01 -------------------------------------------------- >>> Straight-forward. Type = unit diagonal Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 0.40000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 0.29629630E+00 x-vector: 0.35925926E+01 0.00000000E+00 d-vector: 0.75154639E+01 0.33402062E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.24814815E+01 0.37777778E+01 0.50740741E+01 0.63703704E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463954E+02 y-vector: 0.79197001E+01 0.85289078E+01 0.91381155E+01 0.97473232E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.00000000E+00 d-vector: -0.18000000E+01 0.18000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 -0.13333333E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 -0.17888544E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 0.71554175E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.00000000E+00 d-vector: -0.18000000E+01 0.18000000E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 -0.13333333E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 0.53333333E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 -0.17888544E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 0.71554175E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666667E+00 -0.66666667E+00 x-vector: 0.16666667E+01 0.00000000E+00 d-vector: -0.18000000E+01 0.18000000E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: -0.23333333E+01 -0.20000000E+01 -0.16666667E+01 y-vector: 0.43333333E+01 0.46666667E+01 0.50000000E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304952E+01 -0.26832816E+01 -0.22360680E+01 y-vector: 0.58137767E+01 0.62609903E+01 0.67082039E+01 -------------------------------------------------- >>> Attempted row removal but z1>w1. Error. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Input vector is of the form (x1,0). H = I. Type = unit diagonal Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). H = I. Type = unit diagonal Input Data x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,x1). Type = anti-diagonal. Diag elems = 0 Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.30000000E+01 0.00000000E+00 d-vector: 0.90000000E+01 0.40000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: -0.10000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.60000000E+01 y-vector: -0.20000000E+01 -------------------------------------------------- >>> Error: 1/d(2) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(1) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(2) = 0, x(2) = 0. Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(1) = 0, x(1) = 0. Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(2) < 0. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 -0.40000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Scaling test d(1) = d(2) = 1/(1.01*g^2) Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.11347938+308 0.11347938+308 Type -1, Rescaling has taken place Row 1: 0.17272337E-76 -0.17272337E-76 Row 2: 0.17272337E-76 0.17272337E-76 x-vector: 0.11579209E+78 0.00000000E+00 d-vector: 0.14768992-153 0.14768992-153 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.34544674E-76 0.69089348E-76 y-vector: 0.00000000E+00 0.00000000E+00 Output vectors (scaled to match _rot output) x-vector: -0.41981362-153 0.83962724-153 y-vector: 0.00000000E+00 0.00000000E+00 -------------------------------------------------- >>> d(1) = d(2) = huge/four; Type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22250739-307 0.22250739-307 Type 0, H is unit diagonal Off-diagonal elements: 0.10000000E+01 -0.10000000E+01 x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.22471164+308 0.22471164+308 -------------------------------------------------- >>> Scaling test d(1) = d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22471164+308 0.22471164+308 Type -1, Rescaling has taken place Row 1: 0.17272337E-76 -0.17272337E-76 Row 2: 0.17272337E-76 0.17272337E-76 x-vector: 0.11579209E+78 0.00000000E+00 d-vector: 0.74583407-154 0.74583407-154 -------------------------------------------------- >>> Scaling test d(1) = huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22250739-307 0.22471164+308 Type -1, Rescaling has taken place Row 1: 0.10000000E+01 -0.00000000E+00 Row 2: 0.17272337E-76 0.17272337E-76 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.44942328+308 0.14916681-153 -------------------------------------------------- >>> Scaling test d(1) = two*tiny, d(2) = huge/four Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22471164+308 0.22250739-307 Type -1, Rescaling has taken place Row 1: -0.00000000E+00 0.10000000E+01 Row 2: -0.17272337E-76 -0.17272337E-76 x-vector: -0.10000000E+01 0.00000000E+00 d-vector: 0.44942328+308 0.14916681-153 -------------------------------------------------- >>> d(1) = huge/four, d(2) = one; type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.22250739-307 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.10000000E+01 -0.22250739-307 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.44942328+308 0.10000000E+01 -------------------------------------------------- >>> Scaling test d(1) = -huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: -0.22250739-307 0.22471164+308 Type -1, Rescaling has taken place Row 1: 0.10000000E+01 0.00000000E+00 Row 2: 0.17272337E-76 0.17272337E-76 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.44942328+308 0.14916681-153 -------------------------------------------------- >>> d(1) = -two*tiny, d(2) = huge/four, x(1) = 0; Type = anti-diagonal Input Data x-vector: 0.00000000E+00 -0.10000000E+01 d-vector: -0.22471164+308 0.22250739-307 Error condition detected in Rotmg call -------------------------------------------------- >>> d(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal Input Data x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.22471164+308 0.22250739-307 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.44501477-307 0.44942328+308 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res7' then echo shar: will not over-write existing file "'res7'" else cat << "SHAR_EOF" > 'res7' Test Results for generic Fortran 90 version of the modified Givens transformation (single precision) >>> Straight-forward. Type = anti-diagonal. Equal unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.00000000E+00 d-vector: 0.75154643E+01 0.33402061E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 0.91851854E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 0.13333335E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180586E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460391E+01 0.24368310E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Equal non-unit stride Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.00000000E+00 d-vector: 0.75154643E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 0.91851854E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 0.13333335E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 0.25180586E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460391E+01 0.24368310E+01 -------------------------------------------------- >>> Straight-forward. Type = anti-diagonal. Unequal non-unit stride. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.29629630E+00 0.66666669E+00 x-vector: 0.35925927E+01 0.00000000E+00 d-vector: 0.75154643E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.52962961E+01 0.65925927E+01 0.78888888E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 Output vectors (scaled to match _rot output) x-vector: 0.14519450E+02 0.18073162E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460391E+01 -------------------------------------------------- >>> Straight-forward. Type = unit diagonal Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 0.40000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 0.29629630E+00 x-vector: 0.35925927E+01 0.00000000E+00 d-vector: 0.75154643E+01 0.33402061E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.24814816E+01 0.37777777E+01 0.50740738E+01 0.63703704E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: 0.68028197E+01 0.10356531E+02 0.13910242E+02 0.17463955E+02 y-vector: 0.79197001E+01 0.85289078E+01 0.91381149E+01 0.97473221E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.00000000E+00 d-vector: -0.18000001E+01 0.18000001E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 -0.13333335E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 -0.17888547E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.00000000E+00 d-vector: -0.18000001E+01 0.18000001E+01 Output vectors x and y: incx = 2, incy = 2 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 -0.13333335E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 -0.17888547E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 -------------------------------------------------- >>> Row removal. H=0. Type = diagonal. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: -0.66666669E+00 -0.66666669E+00 x-vector: 0.16666665E+01 0.00000000E+00 d-vector: -0.18000001E+01 0.18000001E+01 Output vectors x and y: incx = 2, incy = 3 x-vector: -0.23333335E+01 -0.20000000E+01 -0.16666670E+01 y-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 Output vectors (scaled to match _rot output) x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 -------------------------------------------------- >>> Attempted row removal but z1>w1. Error. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: -0.10000000E+01 0.10000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: -0.66666669E+00 0.66666669E+00 x-vector: 0.16666665E+01 0.00000000E+00 d-vector: 0.18000001E+01 -0.18000001E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.43333335E+01 0.46666665E+01 0.50000000E+01 0.53333330E+01 y-vector: 0.23333335E+01 0.20000000E+01 0.16666670E+01 0.13333335E+01 Output vectors (scaled to match _rot output) x-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 y-vector: 0.31304955E+01 0.26832817E+01 0.22360685E+01 0.17888547E+01 -------------------------------------------------- >>> Input vector is of the form (x1,0). H = I. Type = unit diagonal Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). H = I. Type = unit diagonal Input Data x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.00000000E+00 0.00000000E+00 d-vector: 0.40000000E+01 0.90000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.10000000E+01 y-vector: 0.20000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,x1). Type = anti-diagonal. Diag elems = 0 Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.40000000E+01 0.90000000E+01 Type 1, H is unit anti-diagonal Diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.30000000E+01 0.00000000E+00 d-vector: 0.90000000E+01 0.40000000E+01 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: -0.10000000E+01 Output vectors (scaled to match _rot output) x-vector: 0.60000000E+01 y-vector: -0.20000000E+01 -------------------------------------------------- >>> Error: 1/d(2) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(1) = 0. Input Data x-vector: 0.20000000E+01 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(2) = 0, x(2) = 0. Input Data x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.40000000E+01 0.00000000E+00 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: 1/d(1) = 0, x(1) = 0. Input Data x-vector: 0.00000000E+00 0.30000000E+01 d-vector: 0.00000000E+00 0.90000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Error: d(2) < 0. Input Data x-vector: 0.30000000E+01 0.20000000E+01 d-vector: 0.90000000E+01 -0.40000000E+01 Error condition detected in Rotmg call -------------------------------------------------- >>> Scaling test d(1) = d(2) = 1/(2*g^2) Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.42535296E+38 0.42535296E+38 Type -1, Rescaling has taken place Row 1: 0.46566129E-09 -0.46566129E-09 Row 2: 0.46566129E-09 0.46566129E-09 x-vector: 0.42949673E+10 0.00000000E+00 d-vector: 0.54210109E-19 0.54210109E-19 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.93132257E-09 0.18626451E-08 y-vector: 0.00000000E+00 0.00000000E+00 Output vectors (scaled to match _rot output) x-vector: -0.21684043E-18 0.43368087E-18 y-vector: 0.00000000E+00 0.00000000E+00 -------------------------------------------------- >>> d(1) = d(2) = huge/four; Type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.11754945E-37 0.11754945E-37 Type 0, H is unit diagonal Off-diagonal elements: 0.10000000E+01 -0.10000000E+01 x-vector: 0.20000000E+01 0.00000000E+00 d-vector: 0.42535291E+38 0.42535291E+38 -------------------------------------------------- >>> Scaling test d(1) = d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.42535296E+38 0.42535296E+38 Type -1, Rescaling has taken place Row 1: 0.46566129E-09 -0.46566129E-09 Row 2: 0.46566129E-09 0.46566129E-09 x-vector: 0.42949673E+10 0.00000000E+00 d-vector: 0.54210109E-19 0.54210109E-19 -------------------------------------------------- >>> Scaling test d(1) = huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.11754945E-37 0.42535296E+38 Type -1, Rescaling has taken place Row 1: 0.10000000E+01 -0.00000000E+00 Row 2: 0.46566129E-09 0.46566129E-09 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.85070582E+38 0.10842022E-18 -------------------------------------------------- >>> Scaling test d(1) = two*tiny, d(2) = huge/four Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.42535296E+38 0.11754945E-37 Type -1, Rescaling has taken place Row 1: -0.00000000E+00 0.10000000E+01 Row 2: -0.46566129E-09 -0.46566129E-09 x-vector: -0.10000000E+01 0.00000000E+00 d-vector: 0.85070582E+38 0.10842022E-18 -------------------------------------------------- >>> d(1) = huge/four, d(2) = one; type = anti-diagonal Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: 0.11754945E-37 0.10000000E+01 Type 0, H is unit diagonal Off-diagonal elements: 0.10000000E+01 -0.11754945E-37 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: 0.85070582E+38 0.10000000E+01 -------------------------------------------------- >>> Scaling test d(1) = -huge/four, d(2) = two*tiny Input Data x-vector: 0.10000000E+01 -0.10000000E+01 d-vector: -0.11754945E-37 0.42535296E+38 Type -1, Rescaling has taken place Row 1: 0.10000000E+01 0.00000000E+00 Row 2: 0.46566129E-09 0.46566129E-09 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.85070582E+38 0.10842022E-18 -------------------------------------------------- >>> d(1) = -two*tiny, d(2) = huge/four, x(1) = 0; Type = anti-diagonal Input Data x-vector: 0.00000000E+00 -0.10000000E+01 d-vector: -0.42535296E+38 0.11754945E-37 Type 1, H is unit anti-diagonal Diagonal elements: 0.00000000E+00 -0.00000000E+00 x-vector: -0.10000000E+01 0.00000000E+00 d-vector: 0.85070582E+38 -0.23509887E-37 -------------------------------------------------- >>> d(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal Input Data x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.42535296E+38 0.11754945E-37 Type 0, H is unit diagonal Off-diagonal elements: 0.00000000E+00 0.00000000E+00 x-vector: 0.10000000E+01 0.00000000E+00 d-vector: -0.23509887E-37 0.85070582E+38 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res8' then echo shar: will not over-write existing file "'res8'" else cat << "SHAR_EOF" > 'res8' Test Results for generic Fortran 90 version of the standard Givens transformation (single precision) >>> Straight-forward. Equal unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381150E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.14519449E+02 0.18073160E+02 0.21626873E+02 0.25180584E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460386E+01 0.24368310E+01 -------------------------------------------------- >>> Straight-forward. Equal non-unit stride Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381150E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.14519449E+02 0.18073160E+02 0.21626873E+02 0.25180584E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460386E+01 0.24368310E+01 -------------------------------------------------- >>> Straight-forward. Unequal non-unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381150E+00 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.14519449E+02 0.18073160E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460386E+01 -------------------------------------------------- >>> Straight-forward. Input Data x-vector: 0.90000000E+01 0.40000000E+01 c and s values: 0.91381150E+00 0.40613845E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463953E+02 y-vector: 0.79196997E+01 0.85289078E+01 0.91381149E+01 0.97473221E+01 -------------------------------------------------- >>> d(1) < zero -- drop data. No application step. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535596E+00 0.66666669E+00 -------------------------------------------------- >>> d(1) < zero -- drop data. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535596E+00 0.66666669E+00 -------------------------------------------------- >>> d(1) < zero -- drop data. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535596E+00 0.66666669E+00 -------------------------------------------------- >>> d(1) < zero -- drop data. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535596E+00 0.66666669E+00 -------------------------------------------------- >>> d(1) < zero (z1>w1) -- c=0, s=0. Error condition. Input Data x-vector: 0.20000000E+01 0.30000000E+01 c and s values: 0.00000000E+00 0.00000000E+00 Error detected in d_rot -------------------------------------------------- >>> Input vector is of the form (x1,0). Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). Input Data x-vector: 0.00000000E+00 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 -------------------------------------------------- >>> Input vector is of the form (0,x1). Input Data x-vector: 0.00000000E+00 0.90000000E+01 c and s values: 0.00000000E+00 0.10000000E+01 -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'res9' then echo shar: will not over-write existing file "'res9'" else cat << "SHAR_EOF" > 'res9' Test Results for Fortran 77 compatible version of the standard Givens transformation (single precision) >>> Straight-forward. Equal unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381150E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.14519449E+02 0.18073160E+02 0.21626873E+02 0.25180584E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460386E+01 0.24368310E+01 -------------------------------------------------- >>> Straight-forward. Equal non-unit stride Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381150E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.14519449E+02 0.18073160E+02 0.21626873E+02 0.25180584E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460386E+01 0.24368310E+01 -------------------------------------------------- >>> Straight-forward. Unequal non-unit stride. Input Data x-vector: 0.40000000E+01 0.90000000E+01 c and s values: 0.40613845E+00 0.91381150E+00 Output vectors x and y: incx = 2, incy = 3 x-vector: 0.14519449E+02 0.18073160E+02 0.21626873E+02 y-vector: 0.42644539E+01 0.36552460E+01 0.30460386E+01 -------------------------------------------------- >>> Straight-forward. Input Data x-vector: 0.90000000E+01 0.40000000E+01 c and s values: 0.91381150E+00 0.40613845E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: 0.68028193E+01 0.10356531E+02 0.13910242E+02 0.17463953E+02 y-vector: 0.79196997E+01 0.85289078E+01 0.91381149E+01 0.97473221E+01 -------------------------------------------------- >>> d(1) < zero -- drop data. No application step. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535596E+00 0.66666669E+00 -------------------------------------------------- >>> d(1) < zero -- drop data. Equal unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535596E+00 0.66666669E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 -0.17888546E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 -------------------------------------------------- >>> d(1) < zero -- drop data. Equal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535596E+00 0.66666669E+00 Output vectors x and y: incx = 2, incy = 2 x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 -0.17888546E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 0.71554174E+01 -------------------------------------------------- >>> d(1) < zero -- drop data. Unequal non-unit stride. Input Data x-vector: 0.30000000E+01 0.20000000E+01 c and s values: 0.74535596E+00 0.66666669E+00 Output vectors x and y: incx = 2, incy = 3 x-vector: -0.31304955E+01 -0.26832817E+01 -0.22360685E+01 y-vector: 0.58137770E+01 0.62609901E+01 0.67082043E+01 -------------------------------------------------- >>> d(1) < zero (z1>w1) -- c=0, s=0. Error condition. Input Data x-vector: 0.20000000E+01 0.30000000E+01 c and s values: 0.00000000E+00 0.00000000E+00 Error detected in s_rot -------------------------------------------------- >>> Input vector is of the form (x1,0). Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 Output vectors x and y: incx = 1, incy = 1 x-vector: 0.20000000E+01 y-vector: 0.60000000E+01 -------------------------------------------------- >>> Input vector is of the form (0,0). Input Data x-vector: 0.00000000E+00 0.00000000E+00 c and s values: 0.10000000E+01 0.00000000E+00 -------------------------------------------------- >>> Input vector is of the form (0,x1). Input Data x-vector: 0.00000000E+00 0.90000000E+01 c and s values: 0.00000000E+00 0.10000000E+01 -------------------------------------------------- >>> Error: incx <= 0. Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.00000000E+00 0.00000000E+00 Error detected in s_rot -------------------------------------------------- >>> Error: incy <= 0. Input Data x-vector: 0.40000000E+01 0.00000000E+00 c and s values: 0.00000000E+00 0.00000000E+00 Error detected in s_rot -------------------------------------------------- SHAR_EOF fi # end of overwriting check if test -f 'rot_battery.f90' then echo shar: will not over-write existing file "'rot_battery.f90'" else cat << "SHAR_EOF" > 'rot_battery.f90' PROGRAM rot_test ! Test program for rewritten Givens rotation algorithm. ! This test calls the generic Fortran 90 routine. ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision USE givens_rotations, ONLY : rot ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0E0_wp REAL (wp), PARAMETER :: quarter = 0.25E0_wp REAL (wp), PARAMETER :: three = 3.0E0_wp REAL (wp), PARAMETER :: two = 2.0E0_wp REAL (wp), PARAMETER :: zero = 0.0E0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Local Scalars .. REAL (wp) :: c, s INTEGER :: i, incx, incy, k CHARACTER (80) :: title ! .. ! .. Local Arrays .. REAL (wp) :: d(2), x(2), xv(10), yv(10) ! .. ! .. External Subroutines .. EXTERNAL outres ! .. ! .. Intrinsic Functions .. INTRINSIC sqrt ! .. WRITE (6,'(1x,a,/,a)') & ' Test Results for generic Fortran 90 version of the', & ' standard Givens transformation (double precision)' ! The tests are driven via a data file 10 READ (5,'(a)',end=20) title READ (5,*) x(1:2) READ (5,*) d(1:2) ! if k > 0 then there is update vector data READ (5,*) k IF (k>0) THEN READ (5,*) incx, incy ! incx<1 and incy<1 not sensible -- ignore IF (incx<1 .OR. incy<1) GO TO 10 xv = zero yv = zero READ (5,*) (xv(i),i=1,1+(k-1)*incx,incx) READ (5,*) (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(//">>> ",a)') title ! If d(1) is negative in the data file we require row ! removal. This is flagged to d_rot by a negative value ! of k. IF (d(1) < zero) THEN d(1) = -d(1) k = -k - 1 END IF ! calculate the actual values of w1, z1 from the ! separate d and x values. This allows comparison ! of results between standard and modified methods. x = sqrt(d)*x xv(1:1+(k-1)*incx:incx) = sqrt(d(1))*xv(1:1+(k-1)*incx:incx) yv(1:1+(k-1)*incy:incy) = sqrt(d(2))*yv(1:1+(k-1)*incy:incy) WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) CALL rot(x(1),x(2),k,xv,incx,yv,incy,c,s) CALL outres(c,s,k,xv,yv,incx,incy) GO TO 10 20 END PROGRAM rot_test SUBROUTINE outres(c,s,k,xv,yv,incx,incy) ! Output routine for results ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: zero = 0.0_wp ! .. ! .. Scalar Arguments .. REAL (wp) :: c, s INTEGER :: incx, incy, k ! .. ! .. Array Arguments .. REAL (wp) :: xv(*), yv(*) ! .. ! .. Local Scalars .. INTEGER :: i ! .. ! print up results WRITE (6,'(/)') WRITE (6,'(" c and s values:", 2e16.8)') c, s ! error condition is c=s=0 ! print message and exit IF (c==zero .AND. s==zero) THEN WRITE (6,'(/" Error detected in d_rot"/)') ELSE ! k>0 signifies that data transformation has taken place IF (k>0) THEN WRITE (6,'(/"Output vectors x and y: incx = ",i3,", incy = ",i3)') & incx, incy WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) END IF END IF WRITE (6,'(/"--------------------------------------------------")') END SUBROUTINE outres SHAR_EOF fi # end of overwriting check if test -f 'rotm_battery.f90' then echo shar: will not over-write existing file "'rotm_battery.f90'" else cat << "SHAR_EOF" > 'rotm_battery.f90' PROGRAM rotmtest ! Test program for rewritten modified Givens rotation algorithm. ! The calls are to the generic Fortran 90 version. ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision USE givens_rotations, ONLY : rotm ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0E0_wp REAL (wp), PARAMETER :: quarter = 0.25E0_wp REAL (wp), PARAMETER :: three = 3.0E0_wp REAL (wp), PARAMETER :: two = 2.0E0_wp REAL (wp), PARAMETER :: zero = 0.0E0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Local Scalars .. REAL (wp) :: gamsq INTEGER :: i, incx, incy, k CHARACTER (80) :: title ! .. ! .. Local Arrays .. REAL (wp) :: d(2), param(5), x(2), xv(10), yv(10) ! .. ! .. External Subroutines .. EXTERNAL outres ! .. ! .. Intrinsic Functions .. INTRINSIC huge, min, tiny ! .. WRITE (6,'(1x,a,/,a)') & ' Test Results for generic Fortran 90 version of the', & ' modified Givens transformation (double precision)' ! main loop 10 READ (5,'(a)',end=20) title WRITE (6,'(//">>> ",a)') title READ (5,*) x(1:2) READ (5,*) d(1:2) ! The first tests are driven via a data file ! if k > 0 then there is update vector data READ (5,*) k IF (k>0) THEN READ (5,*) incx, incy xv = zero yv = zero READ (5,*) (xv(i),i=1,1+(k-1)*incx,incx) READ (5,*) (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) ! Generate the reciprocals -- we assume that it is the d ! values that are provided. Anyway in most cases these ! start as one so there is only a risk of confusion ! on exit. IF (d(1)/=zero) THEN d(1) = one/d(1) END IF IF (d(2)/=zero) THEN d(2) = one/d(2) END IF CALL rotm(d(1),d(2),x(1),x(2),k,xv,incx,yv, incy,param) CALL outres(param,x,d,k,xv,yv,incx,incy) GO TO 10 ! tests for rescaling and tests for large and small data values 20 CONTINUE gamsq = min(huge(one),one/tiny(one))*quarter ! scaling test d(1) = d(2) < 1/g^2 x(1) = one x(2) = -one d(1) = one/(two*gamsq) d(2) = d(1) d = one/d k = 2 xv(1) = -one xv(2) = two yv(1) = one yv(2) = -two title = 'Scaling test d(1) = d(2) = 1/(2*g^2)' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = d(2) = huge/four; type = anti-diagonal x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = huge(one)*quarter d = one/d k = -1 title = 'd(1) = d(2) = huge/four; Type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = d(2) = two*tiny x(1) = one x(2) = -one d(1) = two*tiny(one) d(2) = two*tiny(one) d = one/d k = -1 title = 'Scaling test d(1) = d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) =huge/four, d(2) = two*tiny x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = two*tiny(one) d = one/d k = -1 title = 'Scaling test d(1) = huge/four, d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = two*tiny, d(2) = huge/four x(1) = one x(2) = -one d(1) = two*tiny(one) d(2) = huge(one)*quarter d = one/d k = -1 title = 'Scaling test d(1) = two*tiny, d(2) = huge/four' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = huge/four, d(2) = one; type = anti-diagonal x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = one d = one/d k = -1 title = 'd(1) = huge/four, d(2) = one; type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = -huge/four, d(2) = two*tiny x(1) = one x(2) = -one d(1) = -huge(one)*quarter d(2) = two*tiny(one) d = one/d k = -1 title = 'Scaling test d(1) = -huge/four, d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = -two*tiny, d(2) = huge/four, x(1) = 0; type = anti-diagonal x(1) = zero x(2) = -one d(1) = -two*tiny(one) d(2) = huge(one)*quarter d = one/d k = -1 title = & 'd(1) = -two*tiny, d(2) = huge/four, x(1) = 0; Type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = -two*tiny, d(2) = huge/four, x(2) = 0; type = diagonal x(1) = one x(2) = zero d(1) = -two*tiny(one) d(2) = huge(one)*quarter d = one/d k = -1 title = 'd(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) END PROGRAM rotmtest SUBROUTINE outres(param,x,d,k,xv,yv,incx,incy) ! Output routine for the returned values ! .. Use Statements .. USE import_kinds, ONLY : wp => double_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Scalar Arguments .. INTEGER :: incx, incy, k ! .. ! .. Array Arguments .. REAL (wp) :: d(2), param(5), x(2), xv(*), yv(*) ! .. ! .. Local Scalars .. INTEGER :: i ! .. ! .. Intrinsic Functions .. INTRINSIC abs, int, sqrt ! .. ! print up results WRITE (6,'(/)') SELECT CASE (int(param(1))) CASE (rescaled) WRITE (6,'(" Type -1, Rescaling has taken place")') WRITE (6,'(" Row 1: ", 2e16.8)') param(2), param(4) WRITE (6,'(" Row 2: ", 2e16.8/)') param(3), param(5) WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') one/d(1:2) CASE (sltc) WRITE (6,'(" Type 0, H is unit diagonal")') WRITE (6,'(/" Off-diagonal elements:", 2e16.8)') param(3), param(4) WRITE (6,'( " x-vector: ", 2e16.8)') x(1:2) WRITE (6,'( " d-vector: ", 2e16.8)') one/d(1:2) CASE (clts) WRITE (6,'(" Type 1, H is unit anti-diagonal")') WRITE (6,'(/" Diagonal elements: ", 2e16.8)') param(2), param(5) WRITE (6,'( " x-vector: ", 2e16.8)') x(1:2) WRITE (6,'( " d-vector: ", 2e16.8)') one/d(1:2) CASE (error) WRITE (6,'(" Error condition detected in Rotmg call")') WRITE (6,'(/"--------------------------------------------------")') RETURN END SELECT ! k>0 if data has been transformed IF (k>0) THEN WRITE (6,'(/"Output vectors x and y: incx = ",i3,", incy = ",i3)') & incx, incy WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) WRITE (6,'(/"Output vectors (scaled to match _rot output)")') xv(1:1+(k-1)*incx:incx) = xv(1:1+(k-1)*incx:incx)/sqrt(abs(d(1))) yv(1:1+(k-1)*incy:incy) = yv(1:1+(k-1)*incy:incy)/sqrt(abs(d(2))) WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(/"--------------------------------------------------")') END SUBROUTINE outres SHAR_EOF fi # end of overwriting check if test -f 's_rot_battery.f90' then echo shar: will not over-write existing file "'s_rot_battery.f90'" else cat << "SHAR_EOF" > 's_rot_battery.f90' PROGRAM s_rot_test ! Test program for rewritten Givens rotation algorithm. ! This routine combines the two original blas 1 routines ! srot and srotg. ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision USE givens_rotations, ONLY : s_rot ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0E0_wp REAL (wp), PARAMETER :: quarter = 0.25E0_wp REAL (wp), PARAMETER :: three = 3.0E0_wp REAL (wp), PARAMETER :: two = 2.0E0_wp REAL (wp), PARAMETER :: zero = 0.0E0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Local Scalars .. REAL (wp) :: c, s INTEGER :: i, incx, incy, k CHARACTER (80) :: title ! .. ! .. Local Arrays .. REAL (wp) :: d(2), x(2), xv(10), yv(10) ! .. ! .. External Subroutines .. EXTERNAL outres ! .. ! .. Intrinsic Functions .. INTRINSIC sqrt ! .. WRITE (6,'(1x,a,/,a)') & ' Test Results for Fortran 77 compatible version of the', & ' standard Givens transformation (single precision)' ! The tests are driven via a data file 10 READ (5,'(a)',end=20) title WRITE (6,'(//">>> ",a)') title READ (5,*) x(1:2) READ (5,*) d(1:2) ! if k > 0 then there is update vector data READ (5,*) k IF (k>0) THEN READ (5,*) incx, incy xv = zero yv = zero IF (incx>0 .AND. incy>0) THEN READ (5,*) (xv(i),i=1,1+(k-1)*incx,incx) READ (5,*) (yv(i),i=1,1+(k-1)*incy,incy) END IF END IF ! If d(1) is negative in the data file we require row ! removal. This is flagged to s_rot by a negative value ! of k. IF (d(1)0 .AND. incy>0) THEN xv(1:1+(k-1)*incx:incx) = sqrt(d(1))*xv(1:1+(k-1)*incx:incx) yv(1:1+(k-1)*incy:incy) = sqrt(d(2))*yv(1:1+(k-1)*incy:incy) END IF WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) CALL s_rot(x(1),x(2),k,xv,incx,yv,incy,c,s) CALL outres(c,s,k,xv,yv,incx,incy) GO TO 10 20 END PROGRAM s_rot_test SUBROUTINE outres(c,s,k,xv,yv,incx,incy) ! Output routine for results ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: zero = 0.0_wp ! .. ! .. Scalar Arguments .. REAL (wp) :: c, s INTEGER :: incx, incy, k ! .. ! .. Array Arguments .. REAL (wp) :: xv(*), yv(*) ! .. ! .. Local Scalars .. INTEGER :: i ! .. ! .. Intrinsic Functions .. INTRINSIC abs ! .. ! print up results WRITE (6,'(/)') WRITE (6,'(" c and s values:", 2e16.8)') c, s ! error condition is c=s=0 ! print message and exit IF (c==zero .AND. s==zero) THEN WRITE (6,'(/" Error detected in s_rot"/)') ELSE ! Fix k value when row removal is being flagged IF (k<0) THEN k = abs(k) - 1 END IF ! k>0 signifies that data transformation has taken place IF (k>0) THEN WRITE (6,'(/"Output vectors x and y: incx = ",i3,", incy = ",i3)') & incx, incy WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) END IF END IF WRITE (6,'(/"--------------------------------------------------")') END SUBROUTINE outres SHAR_EOF fi # end of overwriting check if test -f 's_rotm_battery.f90' then echo shar: will not over-write existing file "'s_rotm_battery.f90'" else cat << "SHAR_EOF" > 's_rotm_battery.f90' PROGRAM s_rotmtest ! Test program for rewritten modified Givens rotation algorithm. ! This routine combines the two original blas 1 routines ! srotm and srotmg. ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision USE givens_rotations, ONLY : s_rotm ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0E0_wp REAL (wp), PARAMETER :: quarter = 0.25E0_wp REAL (wp), PARAMETER :: three = 3.0E0_wp REAL (wp), PARAMETER :: two = 2.0E0_wp REAL (wp), PARAMETER :: zero = 0.0E0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Local Scalars .. REAL (wp) :: gamsq INTEGER :: i, incx, incy, k CHARACTER (80) :: title ! .. ! .. Local Arrays .. REAL (wp) :: d(2), param(5), x(2), xv(10), yv(10) ! .. ! .. External Subroutines .. EXTERNAL outres ! .. ! .. Intrinsic Functions .. INTRINSIC huge, min, tiny ! .. WRITE (6,'(1x,a,/,a)') & ' Test Results for Fortran 77 compatible version of the', & ' modified Givens transformation (single precision)' ! The first tests are driven via a data file ! if k > 0 then there is update vector data 10 READ (5,'(a)',end=20) title WRITE (6,'(//">>> ",a)') title READ (5,*) x(1:2) READ (5,*) d(1:2) ! if k > 0 then there is update vector data READ (5,*) k IF (k>0) THEN READ (5,*) incx, incy xv = zero yv = zero READ (5,*) (xv(i),i=1,1+(k-1)*incx,incx) READ (5,*) (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) ! Generate the reciprocals -- we assume that it is the d ! values that are provided. Anyway in most cases these ! start as one so there is only a risk of confusion ! on exit. IF (d(1)/=zero) THEN d(1) = one/d(1) END IF IF (d(2)/=zero) THEN d(2) = one/d(2) END IF CALL s_rotm(d(1),d(2),x(1),x(2),k,xv,incx,yv,incy,param) CALL outres(param,x,d,k,xv,yv,incx,incy) GO TO 10 ! tests for rescaling and tests for large and small data values 20 CONTINUE gamsq = min(huge(one),one/tiny(one))*quarter ! scaling test d(1) = d(2) < 1/g^2 x(1) = one x(2) = -one d(1) = one/(two*gamsq) d(2) = d(1) xv(1) = -one xv(2) = two yv(1) = one yv(2) = -two k = 2 d = one/d title = 'Scaling test d(1) = d(2) = 1/(1.01E0_wp*g^2)' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),k,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = d(2) = huge/four; type = anti-diagonal k = -1 x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = huge(one)*quarter d = one/d title = 'd(1) = d(2) = huge/four; Type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = d(2) = two*tiny x(1) = one x(2) = -one d(1) = two*tiny(one) d(2) = two*tiny(one) d = one/d title = 'Scaling test d(1) = d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) =huge/four, d(2) = two*tiny x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = two*tiny(one) d = one/d title = 'Scaling test d(1) = huge/four, d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = two*tiny, d(2) = huge/four x(1) = one x(2) = -one d(1) = two*tiny(one) d(2) = huge(one)*quarter d = one/d title = 'Scaling test d(1) = two*tiny, d(2) = huge/four' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = huge/four, d(2) = one; type = anti-diagonal x(1) = one x(2) = -one d(1) = huge(one)*quarter d(2) = one d = one/d title = 'd(1) = huge/four, d(2) = one; type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! scaling test d(1) = -huge/four, d(2) = two*tiny x(1) = one x(2) = -one d(1) = -huge(one)*quarter d(2) = two*tiny(one) d = one/d title = 'Scaling test d(1) = -huge/four, d(2) = two*tiny' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = -two*tiny, d(2) = huge/four, x(1) = 0; type = anti-diagonal x(1) = zero x(2) = -one d(1) = -two*tiny(one) d(2) = huge(one)*quarter d = one/d title = & 'd(1) = -two*tiny, d(2) = huge/four, x(1) = 0; Type = anti-diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) ! d(1) = -two*tiny, d(2) = huge/four, x(2) = 0; type = diagonal x(1) = one x(2) = zero d(1) = -two*tiny(one) d(2) = huge(one)*quarter d = one/d title = 'd(1) = two*tiny, d(2) = -huge/four, x(2) = 0; Type = diagonal' WRITE (6,'(//">>> ",a)') title WRITE (6,'(/" Input Data")') WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') d(1:2) CALL s_rotm(d(1),d(2),x(1),x(2),-1,xv,1,yv,1,param) CALL outres(param,x,d,k,xv,yv,incx,incy) END PROGRAM s_rotmtest SUBROUTINE outres(param,x,d,k,xv,yv,incx,incy) ! Output routine for the returned values ! .. Use Statements .. USE import_kinds, ONLY : wp => single_precision ! .. ! .. Parameters .. REAL (wp), PARAMETER :: one = 1.0_wp INTEGER, PARAMETER :: clts = 1, error = -1, rescaled = 2, sltc = 0 ! .. ! .. Scalar Arguments .. INTEGER :: incx, incy, k ! .. ! .. Array Arguments .. REAL (wp) :: d(2), param(5), x(2), xv(*), yv(*) ! .. ! .. Local Scalars .. INTEGER :: i ! .. ! .. Intrinsic Functions .. INTRINSIC abs, int, sqrt ! .. ! print up results WRITE (6,'(/)') SELECT CASE (int(param(1))) CASE (rescaled) WRITE (6,'(" Type -1, Rescaling has taken place")') WRITE (6,'(" Row 1: ", 2e16.8)') param(2), param(4) WRITE (6,'(" Row 2: ", 2e16.8/)') param(3), param(5) WRITE (6,'(" x-vector:", 2e16.8)') x(1:2) WRITE (6,'(" d-vector:", 2e16.8)') one/d(1:2) CASE (sltc) WRITE (6,'(" Type 0, H is unit diagonal")') WRITE (6,'(/" Off-diagonal elements:", 2e16.8)') param(3), param(4) WRITE (6,'( " x-vector: ", 2e16.8)') x(1:2) WRITE (6,'( " d-vector: ", 2e16.8)') one/d(1:2) CASE (clts) WRITE (6,'(" Type 1, H is unit anti-diagonal")') WRITE (6,'(/" Diagonal elements: ", 2e16.8)') param(2), param(5) WRITE (6,'( " x-vector: ", 2e16.8)') x(1:2) WRITE (6,'( " d-vector: ", 2e16.8)') one/d(1:2) CASE (error) WRITE (6,'(" Error condition detected in Rotmg call")') WRITE (6,'(/"--------------------------------------------------")') RETURN END SELECT ! k>0 if data has been transformed IF (k>0) THEN WRITE (6,'(/"Output vectors x and y: incx = ",i3,", incy = ",i3)') & incx, incy WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) WRITE (6,'(/"Output vectors (scaled to match _rot output)")') xv(1:1+(k-1)*incx:incx) = xv(1:1+(k-1)*incx:incx)/sqrt(abs(d(1))) yv(1:1+(k-1)*incy:incy) = yv(1:1+(k-1)*incy:incy)/sqrt(abs(d(2))) WRITE (6,'("x-vector:")') WRITE (6,'(5x,4e16.8)') (xv(i),i=1,1+(k-1)*incx,incx) WRITE (6,'("y-vector:")') WRITE (6,'(5x,4e16.8)') (yv(i),i=1,1+(k-1)*incy,incy) END IF WRITE (6,'(/"--------------------------------------------------")') END SUBROUTINE outres SHAR_EOF fi # end of overwriting check if test -f 'sdate_res' then echo shar: will not over-write existing file "'sdate_res'" else cat << "SHAR_EOF" > 'sdate_res' Downdating Test for Givens Transformation 10 random matrices of order 40 are used Errors given in multiples of ulb Average Max Min NEW SG Downdating 0.3423E+03 0.1721E+04 0.2534E+01 NEW MG Downdating 0.3472E+03 0.2066E+04 0.8779E+01 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'Givens_Rotations.f90' then echo shar: will not over-write existing file "'Givens_Rotations.f90'" else cat << "SHAR_EOF" > 'Givens_Rotations.f90' MODULE givens_rotations ! .. Use Statements .. USE import_kinds, ONLY : sp => single_precision, dp => double_precision ! .. ! .. Generic Interface Blocks .. INTERFACE rot MODULE PROCEDURE s_rot MODULE PROCEDURE d_rot END INTERFACE INTERFACE rotm MODULE PROCEDURE s_rotm MODULE PROCEDURE d_rotm END INTERFACE ! .. CONTAINS ! Modified by R. Hanson 24 April 2002. Use more stable downdating. ! Ref. is from Bjorck, Stewart and Stewart, and Chambers. ! !!!!!!!!!!!!!!!!!!!!!!!!! Documentation !!!!!!!!!!!!!!!!!!!!! ! This module contains a set of routines that efficiently implement ! the calculation and application of a Givens plane rotation using ! both the straightforward and modified methods. ! These routines are intended to replace the blas 1 routines ! drot, drotm, drotg, drotmg, srot, srotm, srotg and srotmg. ! Six user callable routines have been provided. ! Calls to d_rot are intended to be direct replacements for a pair of ! calls to the blas level 1 routines drotg and drot; i.e., the ! calculation of the rotation parameters and the application of the ! rotation to a pair of rows has been amalgamated in a single call. ! s_rot replaces calls to srotg and srot. ! In a similar fashion d_rotm and s_rotm replace pairs of calls to drotg ! and drotg, and srotm and srotmg for the modified version of the ! algorithm. ! The routine *_drot and *_drotm mirror the original blas routines in ! their calling sequences. The data to be transformed is passed to the ! routines as a pair of one dimensional, assumed size arrays and the incx ! and incy stride parameters are also available. These routines may be ! used in Fortran 77 codes that use an array element to define the start ! of an array. ! In addition we have provided a pair of generic routines rot and rotm ! that provide the functionality of the *_rot and *_rotm routines ! but allow the same name to be used for either double or single ! precision data. These routines use assumed size arrays to store ! the row data. ! Experiments were conducted using assumed shape arrays designed to ! shorten the calling sequences by allowing the effect of the use of ! strides to be obtained by using a Fortran 90 array slice; for example ! x(1:1+(k-1)*incx:incx). Also the use of assumed shape arrays would ! allow the length of the data matrices to be determined internally using ! the size function. Unfortunately these codes proved to be extrememly ! inefficient at run time and were abandoned. ! The generic routines are implemented as wrappers to the ! *_rot and *_rotm routines. ! In line with the original blas routines both single and double ! precision versions of the routines are available. We have used a ! consistent naming convention for arguments throughout the routines. ! !!!!!!!!!!!!!!!!!!!!!!!!! Routine Descriptions !!!!!!!!!!!!!!!!!!!!!!! ! S_ROT and D_ROT ! --------------- ! construct a Givens transformation, i.e, they compute ! c = cos(theta) and s = sin(theta) such that either ! ( c s) ( w1 ) ( r ) ! ( ) ( ) = ( ) ! (-s c) ( z1 ) ( 0 ) ! or, when dropping data, they compute the hyperbolic ! transformation matrix of the form ! ( c is) ! (-is c) ! where i is the imaginary unit and this unit is implicitly applied. ! When data is dropped the condition, ! w1**2-z1**2 = (w1-z1)*(w1+z1) > 0 ! is required. ! These routines then apply the computed transformation to the data ! ( x_i ) ! ( ) ! ( y_i ) ! where the ith element of the vector x is ! 1 + (i-1)*incx, if incx >= 0, ! and similarly for y and incy. ! Note that these routines do not allow negative values of incx ! and incy, nor do they compute the z value returned by *rotg ! In addition these routines also allow for row removal which ! was only available with the modified Givens routines *ROTMG ! in the original Blas level 1 library. ! S_ROTM and D_ROTM ! ----------------- ! construct the modified Givens rotation matrix H which ! zeros the second component of the 2-vector ! (w1,z1) = (sqrt(d1)*x1, sqrt(d2)*x2)' ! where d1 and d2 are the reciprocals of the scaling factors ! described in the introduction to the paper and appearing ! in the original blas level 1 routines. ! Row removal mathematically requires d2<0 but this is flagged ! in our code by setting d1<0 for efficiency reasons. ! When data is dropped the condition, ! w1**2-z1**2 = (w1-z1)*(w1+z1) > 0 ! is required. ! and apply it to the 2 x N ! ( x' ) ! matrix ( ), where the elements of the x-vector used are ! ( y' ) ! 1 + (i-1)*incx, if incx >= 0, ! and similarly for y and incy. ! If the value of ! NOTE that these routines require the reciprocal of the d values where ! the original blas routines used the d values directly. Also these ! routines do not allow negative values of incx and incy ! ROT and ROTM ! ------------ ! are Fortran 90 generic interfaces to the routines *_ROT and *_ROTM ! respectively. These routines allow the same routine name to be used ! for either single or double precision data. ! In addition, ROTM uses a structure to return the details of the ! transformation matrix used rather than a floating point array. ! !!!!!!!!!!!!!!!!!!!!!!!!! Parameter Descriptions !!!!!!!!!!!!!!!!!!!!!!! ! ! *_ROT, *_ROTM, ROT and ROTM ! --------------------------- ! (In what follows REAL should be replaced by DOUBLE PRECISION ! when calling the D_ROT and D_ROTM routines. Either REAL or DOUBLE ! PRECISION data may be used with the ROT and ROTM routines provided ! that all the real parameters are of the same type in each call) ! W1,Z1 - REAL (ROT only) ! On entry, these define the two vector whose second element is ! to be annihilated. ! On exit, W1 is set to the value of r (see above) and Z1 is ! set to zero. ! K - INTEGER ! On entry, specifies the number of elements (columns) ! in the two data vectors x and y. ! Unchanged on exit. ! X(*),Y(*) - REAL assumed size arrays ! On entry, these contain the row data to be transformed. ! On exit, unless an error is detected in the input data, ! these contain the transformed data. ! INCX - INTEGER ! On entry, INCX specifies the increment parameter used to step ! through the array SX. Must be positive. ! Unchanged on exit. ! INCY - INTEGER ! On entry, INCY specifies the increment parameter used to step ! through the array SY. Must be positive. ! Unchanged on exit. ! C, S - REAL (ROT only) ! On exit, these are elements of the plane rotation or ! hyperbolic matrix. ! X1, X2 - REAL (ROTM only) ! On entry, contains the values that, when combined with the ! scaling factors RD1 and RD2 (below) give the two vector ! (W1, Z1) whose second element is to be annihilated. ! On exit, X1 contains the transformed value and X2 is ! set to zero. ! RD1, RD2 - REAL (ROTM only) ! On entry, contain the scaling values used for the modified ! transformation. These values are the reciprocal of the values ! (D1 and D2) used in the original level 1 blas routines. ! On exit, contain the scaling values for the transformed vector. ! PARAM(5) - REAL array. (ROTM only) ! On exit, the array param is used to control the form of the ! transformation matrix returned by the modified givens ! algorithm by setting param(1) = sflag where ! sflag = -1.0e0 sflag = 0.0e0 sflag = 1.0e0 ! (sh11 sh12) (1.e0 sh12) ( sh11 1.e0) ! H = ( ) ( ) ( ) ! (sh21 sh22) (sh21 1.e0) (-1.e0 sh22) ! In addition if param(1) = -2 then H is set to I. ! Elements 2--5 of param are used to specify the values of ! sh11, sh21, sh12, sh22 respectively. Values of 1.0e0, -1.0e0, ! or 0.0e0, implied by the value of param(1), are stored ! in param but are not used as multipliers in the transformation ! stage. SUBROUTINE s_rotm(rd1,rd2,x1,y1,k,x,incx,y,incy,param) ! Routine to implement the modified Givens transformation ! This routine is equivalent to the two level 1 blas routines ! srotm and srotmg. ! .. Parameters .. REAL (sp), PARAMETER :: one = 1E0_sp REAL (sp), PARAMETER :: quarter = 25E-2_sp REAL (sp), PARAMETER :: zero = 0E0_sp INTEGER, PARAMETER :: clts = 1, error = -1, rescale = 2, slec = 0 ! .. ! .. Scalar Arguments .. REAL (sp), INTENT (INOUT) :: rd1, rd2, x1, y1 INTEGER, INTENT (IN) :: incx, incy, k ! .. ! .. Array Arguments .. REAL (sp), INTENT (OUT) :: param(5) REAL (sp), INTENT (INOUT) :: x(*), y(*) ! .. ! .. Local Scalars .. REAL (sp) :: gam, gamsq, h11, h12, h21, h22, p1, p2, rgam, sx1, sy1, u INTEGER :: i, kx, ky, n ! .. ! .. Intrinsic Functions .. INTRINSIC abs, huge, min, sqrt, tiny ! .. ! Check input parameters; quick return if illegal IF (rd1==zero .OR. rd2<=zero .OR. (k>0 .AND. (incx<=0 .OR. incy<= & 0))) THEN param(1) = error RETURN END IF n = k kx = 1 ky = 1 sx1 = x1 sy1 = y1 h11 = one h21 = zero h12 = zero h22 = one p2 = rd1*sy1 p1 = rd2*sx1 ! This value is dependent on the underlying ! floating-point arithmetic and need be computed once. ! To make the code thread-safe it is done every time. ! When rescaling is needed gam and rgam are computed. gamsq = sqrt(min(huge(one),one/tiny(one))*quarter) ! If |c| >= |s|: IF (p1*sx1>=abs(p2*sy1)) THEN param(1) = slec ! There is no transformation necessary. IF (p2==zero) GO TO 10 h21 = -sy1/sx1 h12 = p2/p1 u = one - h12*h21 sx1 = sx1*u rd1 = rd1*u rd2 = rd2*u ! Check if rescaling is required on either diagonal factor. IF (abs(rd1)+rd2>gamsq) GO TO 20 ! Transformation : (1 h12) ! (h21 1) param(2:5) = (/h11, h21, h12, h22/) ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n p1 = x(i) + param(4)*y(i) y(i) = y(i) + param(3)*x(i) x(i) = p1 END DO ELSE ! Non-unit but equal strides: DO i = 1, n p1 = x(kx) + param(4)*y(kx) y(kx) = y(kx) + param(3)*x(kx) x(kx) = p1 kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n p1 = x(kx) + param(4)*y(ky) y(ky) = y(ky) + param(3)*x(kx) x(kx) = p1 kx = kx + incx ky = ky + incy END DO END IF GO TO 10 ELSE h22 = sx1/sy1 h11 = p1/p2 h12 = one h21 = -one u = one + h11*h22 p1 = rd2*u rd2 = rd1*u sx1 = sy1*u rd1 = p1 param(1) = clts ! Check if rescaling is required on either diagonal factor. IF (abs(rd1)+rd2>gamsq) GO TO 20 ! Transformation : ( h11 1) ! (-1 h22) param(2:5) = (/h11, h21, h12, h22/) ! Perhaps only an interchange and sign change is needed. IF (h11==zero .AND. h22==zero) THEN DO i = 1, n p1 = y(ky) y(ky) = -x(kx) x(kx) = p1 kx = kx + incx ky = ky + incy END DO GO TO 10 END IF ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n p1 = y(i) + param(2)*x(i) y(i) = -x(i) + param(5)*y(i) x(i) = p1 END DO ELSE ! Non-unit but equal strides: DO i = 1, n p1 = y(kx) + param(2)*x(kx) y(kx) = -x(kx) + param(5)*y(kx) x(kx) = p1 kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n p1 = y(ky) + param(2)*x(kx) y(ky) = -x(kx) + param(5)*y(ky) x(kx) = p1 kx = kx + incx ky = ky + incy END DO END IF END IF 10 CONTINUE ! The form of the matrix is defined in param(1) when here. param(2:5) = (/h11, h21, h12, h22/) y1 = zero x1 = sx1 RETURN 20 CONTINUE ! Compute gam and rgam if they may be needed for rescaling. gam = sqrt(gamsq) rgam = one/gam IF (abs(rd1)>gamsq) THEN rd1 = (rd1*rgam)*rgam sx1 = sx1*gam h11 = h11*rgam h12 = h12*rgam param(1) = rescale END IF IF (rd2>gamsq) THEN rd2 = (rd2*rgam)*rgam h21 = h21*rgam h22 = h22*rgam param(1) = rescale END IF ! Apply scaled transformation. This is a rare event, after start. DO i = 1, n p1 = h11*x((i-1)*incx+1) + h12*y((i-1)*incy+1) y((i-1)*incy+1) = h21*x((i-1)*incx+1) + h22*y((i-1)*incy+1) x((i-1)*incx+1) = p1 END DO param(2:5) = (/h11, h21, h12, h22/) GO TO 10 END SUBROUTINE s_rotm SUBROUTINE d_rotm(rd1,rd2,x1,y1,k,x,incx,y,incy,param) ! Routine to implement the modified Givens transformation ! This routine is equivalent to the two level 1 blas routines ! drotm and drotmg. ! .. Parameters .. REAL (dp), PARAMETER :: one = 1E0_dp REAL (dp), PARAMETER :: quarter = 25E-2_dp REAL (dp), PARAMETER :: zero = 0E0_dp INTEGER, PARAMETER :: clts = 1, error = -1, rescale = 2, slec = 0 ! .. ! .. Scalar Arguments .. REAL (dp), INTENT (INOUT) :: rd1, rd2, x1, y1 INTEGER, INTENT (IN) :: incx, incy, k ! .. ! .. Array Arguments .. REAL (dp), INTENT (OUT) :: param(5) REAL (dp), INTENT (INOUT) :: x(*), y(*) ! .. ! .. Local Scalars .. REAL (dp) :: gam, gamsq, h11, h12, h21, h22, p1, p2, rgam, sx1, sy1, u INTEGER :: i, kx, ky, n ! .. ! .. Intrinsic Functions .. INTRINSIC abs, huge, min, sqrt, tiny ! .. ! Check input parameters; quick return if illegal IF (rd1==zero .OR. rd2<=zero .OR. (k>0 .AND. (incx<=0 .OR. incy<= & 0))) THEN param(1) = error RETURN END IF n = k kx = 1 ky = 1 sx1 = x1 sy1 = y1 h11 = one h21 = zero h12 = zero h22 = one p2 = rd1*sy1 p1 = rd2*sx1 ! Set the value of gamsq on first call to the ! routine. This value is dependent on the underlying ! floating-point arithmetic and need be computed once. ! To make the code thread-safe it is done every time. ! When rescaling is needed gam and rgam are computed. gamsq = sqrt(min(huge(one),one/tiny(one))*quarter) ! If |c| >= |s|: IF (p1*sx1>=abs(p2*sy1)) THEN param(1) = slec ! There is no transformation necessary. IF (p2==zero) GO TO 10 h21 = -sy1/sx1 h12 = p2/p1 u = one - h12*h21 sx1 = sx1*u rd1 = rd1*u rd2 = rd2*u ! Check if rescaling is required on either diagonal factor. IF (abs(rd1)+rd2>gamsq) GO TO 20 ! Transformation : (1 h12) ! (h21 1) param(2:5) = (/h11, h21, h12, h22/) ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n p1 = x(i) + param(4)*y(i) y(i) = y(i) + param(3)*x(i) x(i) = p1 END DO ELSE ! Non-unit but equal strides: DO i = 1, n p1 = x(kx) + param(4)*y(kx) y(kx) = y(kx) + param(3)*x(kx) x(kx) = p1 kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n p1 = x(kx) + param(4)*y(ky) y(ky) = y(ky) + param(3)*x(kx) x(kx) = p1 kx = kx + incx ky = ky + incy END DO END IF GO TO 10 ELSE ! exit with an error return if rd1 < 0 at this point IF (rd1gamsq) GO TO 20 ! Transformation : ( h11 1) ! (-1 h22) param(2:5) = (/h11, h21, h12, h22/) ! Perhaps only an interchange and sign change is needed. IF (h11==zero .AND. h22==zero) THEN DO i = 1, n p1 = y(ky) y(ky) = -x(kx) x(kx) = p1 kx = kx + incx ky = ky + incy END DO GO TO 10 END IF ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n p1 = y(i) + param(2)*x(i) y(i) = -x(i) + param(5)*y(i) x(i) = p1 END DO ELSE ! Non-unit but equal strides: DO i = 1, n p1 = y(kx) + param(2)*x(kx) y(kx) = -x(kx) + param(5)*y(kx) x(kx) = p1 kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n p1 = y(ky) + param(2)*x(kx) y(ky) = -x(kx) + param(5)*y(ky) x(kx) = p1 kx = kx + incx ky = ky + incy END DO END IF END IF 10 CONTINUE ! The form of the matrix is defined in param(1) when here. param(2:5) = (/h11, h21, h12, h22/) y1 = zero x1 = sx1 RETURN 20 CONTINUE ! Compute gam and rgam if they may be needed for rescaling. gam = sqrt(gamsq) rgam = one/gam IF (abs(rd1)>gamsq) THEN rd1 = (rd1*rgam)*rgam sx1 = sx1*gam h11 = h11*rgam h12 = h12*rgam param(1) = rescale END IF IF (rd2>gamsq) THEN rd2 = (rd2*rgam)*rgam h21 = h21*rgam h22 = h22*rgam param(1) = rescale END IF ! Apply scaled transformation. This is a rare event, after start. DO i = 1, n p1 = h11*x((i-1)*incx+1) + h12*y((i-1)*incy+1) y((i-1)*incy+1) = h21*x((i-1)*incx+1) + h22*y((i-1)*incy+1) x((i-1)*incx+1) = p1 END DO param(2:5) = (/h11, h21, h12, h22/) GO TO 10 END SUBROUTINE d_rotm SUBROUTINE s_rot(w1,z1,k,x,incx,y,incy,c,s) ! Routine to implement the standard Givens transformation ! This routine is equivalent to the two level 1 blas routines ! srot and srotg. ! .. Parameters .. REAL (sp), PARAMETER :: half = 5E-1_sp REAL (sp), PARAMETER :: one = 1E0_sp REAL (sp), PARAMETER :: quarter = 25E-2_sp REAL (sp), PARAMETER :: zero = 0E0_sp ! .. ! .. Scalar Arguments .. REAL (sp), INTENT (OUT) :: c, s REAL (sp), INTENT (INOUT) :: w1, z1 INTEGER, INTENT (IN) :: incx, incy, k ! .. ! .. Array Arguments .. REAL (sp), INTENT (INOUT) :: x(*), y(*) ! .. ! .. Local Scalars .. REAL (sp) :: a, b, u, v INTEGER :: i, kx, ky, n ! .. ! .. Intrinsic Functions .. INTRINSIC abs, sqrt ! .. ! Test incx and incy are both positive IF (incx<=0 .OR. incy<=0) THEN c = zero s = zero RETURN END IF n = k a = w1 b = z1 ! If transformation is hyperbolic, code will drop data. IF (n<0) GO TO 10 IF (abs(a)>abs(b)) THEN ! Here ABS(w1) > ABS(z1) u = a + a v = b/u ! Note that u and r have the sign of w1 a = sqrt(quarter+v*v) ! Note that c is positive c = half/a s = v*(c+c) a = a*u ELSE ! Here ABS(w1) <= ABS(z1) IF (b/=zero) THEN u = b + b v = a/u ! Note that u and r have the sign of ! z1 (r is immediately stored in a) a = sqrt(quarter+v*v) ! Note that s is positive s = half/a c = v*(s+s) a = a*u ! Here w1 = z1 = 0. ELSE c = one s = zero GO TO 20 END IF END IF ! Possibly B=0, so no arithmetic is needed. IF (c==one .AND. s==zero) GO TO 20 kx = 1 ky = 1 ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n u = c*x(i) + s*y(i) y(i) = -s*x(i) + c*y(i) x(i) = u END DO ELSE ! Non-unit but equal strides: DO i = 1, n u = c*x(kx) + s*y(kx) y(kx) = -s*x(kx) + c*y(kx) x(kx) = u kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n u = c*x(kx) + s*y(ky) y(ky) = -s*x(kx) + c*y(ky) x(kx) = u kx = kx + incx ky = ky + incy END DO END IF GO TO 20 10 CONTINUE n = -(n+1) ! This is an error condition. IF(ABS(a) <= ABS(b)) THEN c = zero s = zero a = zero RETURN ENDIF ! Will have abs(a) > abs(b) here: s=b/a c=(one-s)*(one+s) ! May have c <= 0, even though mathmatically it should be > 0. ! This is also an error condition. IF(c <= zero) THEN c = zero s = zero a = zero RETURN ENDIF c=sqrt(c) a=c*a ! Will need reciprocal of c to apply transformation: b=one/c kx = 1 ky = 1 ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n u = b*(x(i) - s*y(i)) y(i) = -s*u + c*y(i) x(i) = u END DO ELSE ! Non-unit but equal strides: DO i = 1, n u = b*(x(kx) - s*y(kx)) y(kx) = -s*u + c*y(kx) x(kx) = u kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n u = b*(x(kx) - s*y(ky)) y(ky) = -s*u + c*y(ky) x(kx) = u kx = kx + incx ky = ky + incy END DO END IF 20 CONTINUE z1 = zero w1 = a END SUBROUTINE s_rot SUBROUTINE d_rot(w1,z1,k,x,incx,y,incy,c,s) ! Routine to implement the standard Givens transformation ! This routine is equivalent to the two level 1 blas routines ! drot and drotg. ! .. Parameters .. REAL (dp), PARAMETER :: half = 5E-1_dp REAL (dp), PARAMETER :: one = 1E0_dp REAL (dp), PARAMETER :: quarter = 25E-2_dp REAL (dp), PARAMETER :: zero = 0E0_dp ! .. ! .. Scalar Arguments .. REAL (dp), INTENT (OUT) :: c, s REAL (dp), INTENT (INOUT) :: w1, z1 INTEGER, INTENT (IN) :: incx, incy, k ! .. ! .. Array Arguments .. REAL (dp), INTENT (INOUT) :: x(*), y(*) ! .. ! .. Local Scalars .. REAL (dp) :: a, b, u, v INTEGER :: i, kx, ky, n ! .. ! .. Intrinsic Functions .. INTRINSIC abs, sqrt ! .. ! Test incx and incy are both positive IF (incx<=0 .OR. incy<=0) THEN c = zero s = zero RETURN END IF n = k a = w1 b = z1 ! If transformation is hyperbolic, code will drop data. IF (n<0) GO TO 10 IF (abs(a)>abs(b)) THEN ! Here ABS(w1) > ABS(z1) u = a + a v = b/u ! Note that u and r have the sign of W1 a = sqrt(quarter+v*v) ! Note that c is positive c = half/a s = v*(c+c) a = a*u ELSE ! Here ABS(w1) <= ABS(z1) IF (b/=zero) THEN u = b + b v = a/u ! Note that u and r have the sign of ! z1 (r is immediately stored in a) a = sqrt(quarter+v*v) ! Note that s is positive s = half/a c = v*(s+s) a = a*u ! Here W1 = Z1 = 0. ELSE c = one s = zero GO TO 20 END IF END IF ! Possibly B=0, so no arithmetic is needed. IF (c==one .AND. s==zero) GO TO 20 kx = 1 ky = 1 ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n u = c*x(i) + s*y(i) y(i) = -s*x(i) + c*y(i) x(i) = u END DO ELSE ! Non-unit but equal strides: DO i = 1, n u = c*x(kx) + s*y(kx) y(kx) = -s*x(kx) + c*y(kx) x(kx) = u kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n u = c*x(kx) + s*y(ky) y(ky) = -s*x(kx) + c*y(ky) x(kx) = u kx = kx + incx ky = ky + incy END DO END IF GO TO 20 10 CONTINUE n = -(n+1) IF(ABS(a) <= ABS(b)) THEN c = zero s = zero a = zero RETURN ENDIF ! Will have abs(a) > abs(b) here: s=b/a c=(one-s)*(one+s) ! May have c <= 0, even though mathmatically it should be > 0. ! This is also an error condition. IF(c <= zero) THEN c = zero s = zero a = zero RETURN ENDIF c=sqrt(c) a=c*a ! Will need reciprocal of c to apply transformation: b=one/c kx = 1 ky = 1 ! Equal increments are likely: IF (incx==incy) THEN ! Unit strides are usually best: IF (incx==1) THEN DO i = 1, n u = b*(x(i) - s*y(i)) y(i) = -s*u + c*y(i) x(i) = u END DO ELSE ! Non-unit but equal strides: DO i = 1, n u = b*(x(kx) - s*y(kx)) y(kx) = -s*u + c*y(kx) x(kx) = u kx = kx + incx END DO END IF ELSE ! Non-equal strides: DO i = 1, n u = b*(x(kx) - s*y(ky)) y(ky) = -s*u + c*y(ky) x(kx) = u kx = kx + incx ky = ky + incy END DO END IF 20 CONTINUE z1 = zero w1 = a END SUBROUTINE d_rot END MODULE givens_rotations SHAR_EOF fi # end of overwriting check if test -f 'Import_Kinds.f90' then echo shar: will not over-write existing file "'Import_Kinds.f90'" else cat << "SHAR_EOF" > 'Import_Kinds.f90' MODULE import_kinds ! .. Parameters .. INTEGER, PARAMETER :: double_precision = kind(1.0D0) INTEGER, PARAMETER :: single_precision = kind(1.0E0) ! .. ! .. Intrinsic Functions .. INTRINSIC kind ! .. END MODULE import_kinds SHAR_EOF fi # end of overwriting check cd .. cd .. # End of shell archive exit 0