1 / 40

2.1 BOOKKEEPING OPERATIONS

2.1 BOOKKEEPING OPERATIONS. SCOPY and SSWAP. SCOPY call SCOPY (n, x, incx, y, incy) ex. call SCOPY( 15 , x( 10 ),1,y(20),1) y(19+1*i)  x( 9+1*i ) i= 1:15 SSWAP call SSWAP (n, x, incx, y, incy) ex. call SSWAP(10, x(20), 3 ,y(2), 2 )

geri
Télécharger la présentation

2.1 BOOKKEEPING OPERATIONS

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 2.1 BOOKKEEPING OPERATIONS

  2. SCOPY and SSWAP • SCOPY call SCOPY (n, x, incx, y, incy) ex. call SCOPY(15, x(10),1,y(20),1) y(19+1*i)  x(9+1*i) i= 1:15 • SSWAP call SSWAP (n, x, incx, y, incy) ex. call SSWAP(10, x(20),3,y(2),2) x(17+3*i) <--> y(2*i) i= 1:20 • n : the loop from 1 to n x(index) : a vector starts from index incx,incy : increment number

  3. Example 2.1-1 MATRIX COPY subroutine MATCOP ( m, n, B, bdim, A, adim ) integer m, n, bdim, adim real B(bdim,*), A(adim,*) integer k do 10 k = 1,n call SCOPY( m, B(1,k), 1, A(1,k), 1) 10 continue return end

  4. Example 2.1-2 MATRIX TRANSPOSE subroutine TRANSP( n, A, adim ) integer n, adim real A(admin,*) integer k do 10 k= 1, n-1 call SSWAP( n-k, A(k+1,k), 1, A(k,k+1), adim) 10 continue return end

  5. Example 2.1-3 CIRCULANT MATRIX subroutine CIRCUL( n, C, cdim, v ) integer n, cdim real C(cdim,*), v(*) integer k call SCOPY (n,v(1),1,C(1,1),1) do 10 k = 2,n call SCOPY( n-1, C(2,k-1), 1, C(1,k), 1) C(n,k) = C(1,k-1) 10 continue return end

  6. 2.2 VECTOR OPERATIONS

  7. Scaling Vectors • SSCAL To product of the form a*x, where a is a scalar and x is a vector: call SSCAL (n, a, x, incx ) ex. call SSCAL (50, a, x(2),2) x(2i)  a*x(2i) i= 1:50

  8. Adding a Multiple of One Vector to Another • SAXPY Perform the operation y a*x + y,where a is a scalar and x and y are vectors: call SAXPY (n, a, x, incx, y, incy) ex. call SAXPY (5, a, x(3), 2, y(10), 4) y(6+4i)  a*x(1+2i) + y(6+4i) i= 1:5

  9. Inner Products • SDOT To compute the value of inner product: w = SDOT (n, a, x, incx, y, incy) ex. W = call SDOT (50, x, 1, y, 1) w = x(1)y(1) + … + x(50)y(50) ex. W = call SDOT (50, x, 2, y, 1) w = x(1)y(1) + x(3)y(2) + … + x(50)y(50) +2

  10. Example 2.2-1 Univariate least squares subroutine LS( m, a, b, x ) integer m real a(*), b(*), x real SDOT x = SDOT( m, a, 1, a, 1 ) if ( x .NE. 0.0) then x = SDOT( m, a, 1, a, 1 ) call SAXPY( m, -x, a, 1, b, 1) call SSCAL( m, x, a, 1) endif return end

  11. Example 2.2-2 Skew and symmetric Parts subroutine SYMSKU( n, A, adim ) integer n,adim real A(adim, *) real half integer k,kpl,nmk parameter( half = .50e0) do 10 k = 1, n-1 kp1 = k + 1 nmk = n – k call SAXPY( nmk, 1., A(kp1,k), 1, A(k,kp1), adim ) call SSCAL( nmk, half, A(k,kpl),adim ) call SAXPY( nmk, -1., A(k,kp1), adim, A(kp1,k), 1 ) 10 continue return end

  12. Example 2.2-3 Matrix-Vector Multiplication subroutine MMULT( m, n, A, adim, x, y, job ) integer m, n, adim, job real A(adim,*), x(*), y(*) real SDOT integer i,j if (job .EQ. 0) then do 10 i = 1,m y(i) = 0.0 10 continue do 20 j = 1,n call SAXPY( m, x(j), A(1,j), 1, y(1), 1) 20 continue else do 30 j = 1,n y(j) = SDOT ( m, A(1,j), 1, x(1), 1) 30 continue endif return end

  13. Example 2.2-4 Krylov Matrix subroutine KRYLOV( n, A, adim, K, kdim, v, j ) integer n, adim, kdim, j real A(adim,*), K(kdim,*), v(*) integer p call SCOPY( n, v(1), 1, K(1,1), 1 ) do 10 p = 2,j call MMULT( n, n, A, adim, K(1,p-1), k(1,p), 0 ) 10 continue return end

  14. 2.3 NORM COMPUTATIONS

  15. 1-Norms and 2-Norms • SNRM2 & SASUM The 2-norm of a real n-vector x : sw = SNRM2 (n, x, incx) The 1-norm of a real n-vector x : sw = SASUM (n, a, x, incx) ex. sw = SNRM2 (20, x, 1) w  sqrt( |x(1)|2 + … + |x(20)|2 ) ex. sw = SASUM (20, x, 1) w  |x(1)| + … + |x(20)|

  16. Maximal Entry • ISAMAX Find the index of element having max. absolute value: imax = ISAMAX(n, x, incx) ex. If x = (2, -1, 0, 6, -7) imax = ISAMAX(5, x, 1) imax = 5 imax = ISAMAX(2, x(2), 2) imax = 2

  17. Example 2.3-1 Vector Infinity-Norm real function NORMI(x,n) integer n real x(*) integer ISAMAX real abs NORMI = abs( x( ISMAX( n, x, 1) ) ) return end

  18. Example 2.3-2 Frobenius Matrix Norm real function SNORMF( m, n, A, adim ) integer m, n, adim real A(adim,*) real v(2), SNRM2 integer k SNORMF = 0.0 do 10 k = 1,n v(1) = SNORMF v(2) = SNRM2( m, A(1,k), 1 ) SNORMF =SNRM2 ( 2, v(1), 1 ) 10 continue return end

  19. Example 2.3-3 Matrix 1-Norm real function SNORM1( m, n, A, adim ) integer m, n, adim real A(adim,*) real s, SASUM integer k SNORM1 = 0.0 do 10 k = 1,n s = SASUM( m, A(1,k), 1 ) SNORM1 = max( s, SNORM1) 10 continue return end

  20. Example 2.3-4 Householder Vector subroutine HOUSE( n, x ) integer n real x(*) real alfa,beta,x1,SNRM2, abs if (alfa .NE. 0.0) then x1 = x(1) x(1) = x1 + sign(alfa,x1) beta = sqrt( 2.0 * alfa * (alfa + abs(x1)) ) call SSCAL( n, 1.0/beta, x(1), 1 ) else x(1) = 1 endif return end

  21. 2.4 GIVENS ROTATIONS

  22. Constructing a Givens Rotation I • SROTG It can be used to construct a Givens rotation for the purpose of zeroing the second component of a prescribed 2-vector. call SROTG ( a, b, c, s ) • a, c ,and s : • b : if (a>b) then b=s esleif (b>a and c<>0) b=1/c else b=1

  23. Constructing a Givens Rotation II • Consider an example: x= ( 1, 4, 7, 3 ) We have : call SROTG( x(2), x(4), c, s)

  24. Reconstructing a Givens Rotation subroutine ZTOCS( z, c, s ) real z,c,s if ( abs(z) .LT. 1 ) then s = z c = sqrt(1. - s*s) elseif ( abs(z) .GT. 1 ) then c = 1./z s = sqrt(1. - c*c) else c = 0. s = 1. endif return end

  25. Applying a Givens Rotation • SROT call SROT(n, x, incx, y, incy, c, s) means • call SROT(n, A(p,1), adim, A(q,1), adim, c, s) • call SROT(m, A(1,p), 1, A(1,q), 1, c, s)

  26. Example 2.4-1 Zeroing an n-Vector subroutine ZERO(n,x) integer n real x(*) real c,s integer k do 10 k = n-1,1,-1 call SROTG(x(k),x(k+1),c,s) 10 continue return end

  27. Example 2.4-2 Factored Q Times Vector subroutine APPLY( n, x, y ) integer n real x(*), y(*) real c, s integer k do 10 k = n-1, 1, -1 call ZTOCS( x(k+1), c, s ) call SROT( 1, y(k), 1, y(k+1), 1, c, s ) 10 continue return end

  28. Example 2.4-3 Hessenberg QR Factorization subroutine HESS( n, A, adim) integer n, adim real A(adim,*) real c, s integer j, jp1 do 20 j = 1, n-1 jp1 = j + 1 call SROTG( A(i,j), A(jp1, j), c, s ) call SROT( n-j, A(j,jp1), adim, & A(jp1,jp1), adim, c, s) 20 continue return end

  29. 2.5 DOUBLE PRECISION AND COMPLEX VERSIONS

  30. Copying a Vector • call SCOPY(n, sx, incx, sy, incy) call DCOPY(n, dx, incx, dy, incy) call CCOPY(n, cx, incx, cy, incy) call ZCOPY(n, zx, incx, zy, incy) • s : real single precision d : real double precision c : complex single precision z : double precision

  31. Interchanging Vectors • call SSWAP(n, sx, incx, sy, incy) call DSWAP(n, dx, incx, dy, incy) call CSWAP(n, cx, incx, cy, incy) call ZSWAP(n, zx, incx, zy, incy) • s : real single precision d : real double precision c : complex single precision z : double precision

  32. Multiplying a Vector by a Scalar • call SSCAL(n, sa, sx, incx) call DSCAL(n, da, dx, incx) call CSCAL(n, ca, cx, incx) call ZSCAL(n, za, zx, incx) • Complex case call CSSCAL(n, sa, cx, incx) call ZDSCAL(n, da, zx, incx)

  33. Inner Product • sw = SDOT(n, sx, incx, sy, incy) dw = DDOT(n, dx, incx, dy, incy) • complex case for xHy : cw = CDOTC(n, cx, incx, cy, incy) zw = ZDOTC(n, zx, incx, zy, incy) for xTy : cw = CDOTU(n, cx, incx, cy, incy) zw = ZDOTU(n, zx, incx, zy, incy)

  34. Adding a Multiple of One Vector to Another • call SAXPY(n, sa, sx, incx, sy, incy) call DAXPY(n, da, dx, incx, dy, incy) call CAXPY(n, ca, cx, incx, cy, incy) call ZAXPY(n, za, zx, incx, zy, incy)

  35. 2-Norm • sw = SNRM2(n, sx, incx) dw = DNRM2(n, dx, incx) sw = SCNRM2(n, cx, incx) dw = DZNRM2(n, zx, incx)

  36. 1-Norm • sw = SASUM(n, sx, incx) dw = DASUM(n, dx, incx) sw = SCASUM(n, cx, incx) dw = DZASUM(n, zx, incx)

  37. Maximal Entry • imax = ISAMAX(n, sx, incx) imax = IDAMAX(n, dx, incx) imax = ICAMAX(n, cx, incx) imax = IZAMAX(n, zx, incx)

  38. Constructing a Givens Rotation • call SROTG( sa, sb, sc, ss) call DROTG( da, db, dc, ds)

  39. Applying a Givens Rotation • call SROT(n, sx, incx, sy, incy, sc, ss) call DROT(n, dx, incx, dy, incy, dc, ds) call CSROT(n, cx, incx, cy, incy, sc, ss) call ZDROT(n, zx, incx, zy, incy, dc, ds)

  40. Example 2.5-1 Complex Rotation subroutine SCROTG( ca, cb, sc1, ss1, sc2, ss2) complex ca,cb real sc1,ss1,sc2,ss2 real a1, a2, b1, b2 real real aimag a1 = real(ca) a2 = aimag(ca) b1 = real(cb) b2 = aimag(cb) call SROTG( a1, a2, sca, ssa ) call SROTG( b1, b2, scb, ssb ) call SROTG( a1, b1, sc1, ss1 ) sc2 = sca*scb + ssa*ssb ss2 = ssa*scb - sca*ssb return end

More Related