! ********************************************************************* ! * * ! * subroutine dcvit3 * ! * * ! ********************************************************************* ! Double Precision Version 2.0 ! Written by Gordon A. Fenton, TUNS, Aug. 2, 1993 ! ! PURPOSE computes the initial k1*k2*k3 x k1*k2*k3 covariance matrix (and a ! 27x27 submatrix) between local averages. Used by LAS3G. ! ! ! This routine forms the covariance matrix R0 between local averages arranged ! on a k1 x k2 x k3 grid using a user provided variance function "dvarfn" that ! returns the variance of a local average (note that this is the point ! variance sigma^2 times the traditionally "variance" function, as defined ! by Vanmarcke in "Random Fields", MIT Press, 1984). The grid and element ! numbering schemes appear as follows (for k1 = k2 = 5, k3 = 1 plane shown) ! ! | Dx | ! --- -------------------------- ! Dy | 21 | 22 | 23 | 24 | 25 | ! --- -------------------------- ! | 16 | 17 | 18 | 19 | 20 | ! -------------------------- ---------------- ! | 11 | 12 | 13 | 14 | 15 | | 7 | 8 | 9 | ! -------------------------- ---------------- ! | 6 | 7 | 8 | 9 | 10 | | 4 | 5 | 6 | ! -------------------------- ---------------- ! | 1 | 2 | 3 | 4 | 5 | | 1 | 2 | 3 | ! -------------------------- ---------------- ! ! R0 Array R Array ! ! with the k1 x k2 array numbering on the left and a basic 3 x 3 subarray ! shown on the right. In the depth direction (not shown), the numbering ! continues; ie the element immediately behind element 1 on the left is ! element number 26, etc. If we call Z_i the local average of a random process ! over the i'th cell, then for E[Z] = 0, the elements of R0 are defined by ! ! R0(i,j) = E[Z_i*Z_j] ! ! Note that the elements of R are simply extracted directly from the ! appropriate elements of R0 (for example, R(1,5) = R0(1,7)) if k1, k2 and k3 ! are greater than 2. Note also that since both R0 and R are symmetric, ONLY ! THE UPPER TRIANGULAR VALUES are stored herein. Finally note that ! the random process is assumed to be quadrant symmetric (that is the ! covariance function is even in each individual lag component) to reduce ! the total number of covariance calculations. ! ! Arguments to this routine are as follows; ! ! dvarfn external double precision function which returns the ! variance of a local average of size X x Y x Z. This function ! is referenced using the call "V = dvarfn(X,Y,Z)" (other ! parameters to the function must be passed via common blocks). ! ! R0 double precision array of size at least k1*k2*k3 x k1*k2*k3 ! which on output will contain the covariance matrix between ! the local average elements shown above on the left. Note ! that since R0 is symmetric, only the upper triangular values ! are actually stored (compatible with DCHOL2). (output) ! ! iq leading dimension of the array R0 as specified in the calling ! routine. (input) ! ! R double precision array of size at least 27 x 27 which on output ! will contain the covariance matrix between the local average ! elements shown above on the right. Note that since R is ! symmetric, only the upper triangular values are actually ! stored (compatible with DSIFA). (output) ! ! ir leading dimension of the array R as specified in the calling ! routine. (input) ! ! k1 number of local average cells in the x direction ! (horizontally). (input) ! ! k2 number of local average cells in the y direction ! (vertically). (input) ! ! k3 number of local average cells in the z direction ! (depth). (input) ! ! Dx x-dimension of the local average cells. (input) ! ! Dy y-dimension of the local average cells. (input) ! ! Dz z-dimension of the local average cells. (input) ! ! Required: ! 1) from libGAFsim: DCVAA3 !---------------------------------------------------------------------------- subroutine dcvit3( dvarfn, R0, iq, R, ir, k1, k2, k3, Dx, Dy, Dz ) implicit real*8 (a-h,o-z) dimension R0(iq,*), R(ir,*) external dvarfn ! first form the essential elements of R0 ii = 0 do 30 k = 1, k3 tz = dble(k-1) do 20 j = 1, k2 ty = dble(j-1) do 10 i = 1, k1 ii = ii + 1 tx = dble(i-1) R0(1,ii) = dcvaa3( dvarfn, Dx, Dy, Dz, tx, ty, tz ) 10 continue 20 continue 30 continue ! now distribute into the upper triangle kk = k1*k2*k3 k12 = k1*k2 do 50 n = 2, kk n3 = (n-1)/k12 n2 = (n - k12*n3 - 1)/k1 do 40 m = 2, n m3 = (m-1)/k12 m2 = (m - k12*m3 - 1)/k1 i3 = m3 - n3 i2 = m2 - n2 i1 = (m-n) - k1*i2 - k12*i3 k = 1 + iabs(i1) + k1*iabs(i2) + k12*iabs(i3) R0(m,n) = R0(1,k) 40 continue 50 continue ! form R (27 x 27) if( k1 .lt. 3 .or. k2 .lt. 3 .or. k3 .lt. 3 ) then ! find essential elements of R ii = 0 do 80 k = 1, 3 tz = dble(k-1) do 70 j = 1, 3 ty = dble(j-1) do 60 i = 1, 3 tx = dble(i-1) ii = ii + 1 R(1,ii) = dcvaa3( dvarfn, Dx, Dy, Dz, tx, ty, tz ) 60 continue 70 continue 80 continue ! and distribute into upper tri. do 100 n = 2, 27 n3 = (n-1)/9 n2 = (n - 9*n3 - 1)/3 do 90 m = 2, n m3 = (m-1)/9 m2 = (m - 9*m3 - 1)/3 i3 = m3 - n3 i2 = m2 - n2 i1 = (m-n) - 3*i2 - 9*i3 k = 1 + iabs(i1) + 3*iabs(i2) + 9*iabs(i3) R(m,n) = R(1,k) 90 continue 100 continue else ! or extract R from R0 do 120 n = 1, 27 n3 = (n-1)/9 n2 = (n - 9*n3 - 1)/3 do 110 m = 1, n m3 = (m-1)/9 m2 = (m - 9*m3 - 1)/3 i3 = m3 - n3 i2 = m2 - n2 i1 = (m-n) - 3*i2 - 9*i3 k = 1 + iabs(i1) + k1*iabs(i2) + k12*iabs(i3) R(m,n) = R0(1,k) 110 continue 120 continue endif return end