f_diffusion.f Wed Oct 16 07:59:49 2002 page 1 FILE "f_diffusion.f" 1 program f_diffusion 2 parameter (ndim =2048, sideleft = 10.0, sideright = 40.0, 3 * sidetop = 20.0, sidebottom = 20.0) 4 c 5 c the compile option needed to bring in the LINPACK library are: 6 c f90 -o f_diffusion f_diffusion.f -xlic_lib=sunperf 7 c if you add additional options at the compile line, be sure 8 c that the -xlic_lib option is the final one 9 c 10 c or 11 c f90 -o f_diffusion f_diffusion.f sgeco.o sgesl.o bla s.o 12 c 13 c 14 c set up a linear system A X = b with a maximum dimension on A 15 c given by the parameter ndim. The vector B is called the Right 16 c Hand Side (RHS, for short), the vector X is called the solution. 17 c 18 c Fortran is NOT case sensitive, so the Capital Letters are used 19 c only for emphasis 20 c 21 22 real A(ndim,ndim) 23 c 24 c X stores the computed solution; this storage also holds B initially 25 real X(ndim) 26 c 27 c work is the required floating point workspace 28 real work(ndim) 29 c 30 c ipvt is the required integer workspace 31 integer ipvt(ndim), row, col 32 c -------------------------------------------------------------------- 33 c Two dimensional heat diffusion problem 34 c -------------------------------------------------------------------- 35 c2345678901234567890 36 99999 continue 37 print*,'enter number of rows, K ' 38 print*,' suggested values are : 10, 20, 40 ...' 39 read*,K 40 print*,'enter number of columns, M' 41 print*,' suggested values are : 10, 20, 40 ...' 42 read*,M 43 N = K * M 44 if (N .gt. ndim) then 45 print*, 46 * ' enter values of K and M so their product < ',ndim 47 go to 99999 48 endif 49 50 c 51 c ------------------------------------- 52 c initialize entire matrix to zero using array sections 53 c 54 A(1:N,1:N) = 0.0 55 c 56 c set up the N by N block-structured matrix matrix 57 c we store in a full matrix for experimentation using the Linpack 58 c software SGECO, SGESL instead of the more efficient storage as 59 c just the nono-zero entries (whose notation is slightly complicated) 60 c 61 c for row=1, we have a very special "boundary case" 62 c234567890 63 index = 1 64 A(index,index) = 4.0 f_diffusion.f Wed Oct 16 07:59:49 2002 page 2 65 A(index,index+1) = -1.0 66 A(index,index+M) = -1.0 67 do i=2,M-1 68 index = index + 1 69 A(index,index) = 4.0 70 A(index,index-1) = -1.0 71 A(index,index+1) = -1.0 72 A(index,index+M) = -1.0 73 end do 74 index = index + 1 75 A(index,index) = 4.0 76 A(index,index-1) = -1.0 77 A(index,index+M) = -1.0 78 c 79 c for row = 2 ... K-1, we have the standard pattern 80 c 81 c for col = 2 ... m-1, we have the standard pattern 82 c 83 do row = 2,k-1 84 index = index+1 85 c first grid point of each row is special 86 A(index,index) = 4.0 87 A(index,index+1) = -1.0 88 A(index,index-M) = -1.0 89 A(index,index+M) = -1.0 90 do col = 2,m-1 91 index = index+1 92 A(index,index) = 4.0 93 A(index,index-1) = -1.0 94 A(index,index+1) = -1.0 95 A(index,index-M) = -1.0 96 A(index,index+M) = -1.0 97 end do 98 c last grid point of each row is special 99 index = index+1 100 A(index,index) = 4.0 101 A(index,index-1) = -1.0 102 A(index,index+M) = -1.0 103 A(index,index-M) = -1.0 104 end do 105 c 106 c last row is very special boundary case 107 index = index+1 108 A(index,index) = 4.0 109 A(index,index-M) = -1.0 110 A(index,index+1) = -1.0 111 do col=2,M-1 112 index = index+1 113 A(index,index) = 4.0 114 A(index,index-1) = -1.0 115 A(index,index+1) = -1.0 116 A(index,index-M) = -1.0 117 end do 118 index = index+1 119 A(index,index) = 4.0 120 A(index,index-1) = -1.0 121 A(index,index-M) = -1.0 122 c 123 c ------------------------------------ 124 c set up the RHS vector - mostly zeros 125 X(1:N) = 0.0 126 c boundary conditions 10 + 20 127 X(1) = sideleft + sidebottom 128 do index = 2,m-1 129 c boundary condition on bottom row = 20 degrees 130 X(index) = sidebottom f_diffusion.f Wed Oct 16 07:59:49 2002 page 3 131 end do 132 c boundary conditions 20 + 40 degrees 133 X(M) = sideright + sidetop 134 index = M 135 do row= 2,k-1 136 index = index + 1 137 c boundary condition on the right side = 40 degrees 138 X(index) = sideleft 139 index = index+M-1; 140 c boundary condition on left side = 10 degrees 141 X(index) = sideright 142 end do 143 c very top row 144 index = index + 1 145 X(index) = sideleft + sidetop; 146 do col = 2,m-1 147 index = index + 1 148 c boundary condition on the top = 20 degrees 149 X(index) = sidetop 150 end do 151 c boundary conditions 40 + 20 152 index = index + 1; 153 X(index) = sidetop + sideright 154 c 155 c the software is efficient in its use of memory - and will take as 156 c input the vector of the Right Hand Hand, and overwrite with the 157 c solution vector, X. 158 c ----------------------------------------- 159 c compute the lu decomposition of the matrix A using sgeco from the 160 c LINPACK collection of stable, efficient mathematical software 161 c 162 c the matrix A is a by nm1 by nm1 and its "declared dimension", 163 c which allocates the actual storage in memory is given by ndim 164 c 165 call sgeco (A,ndim,N,ipvt,rcond,work) 166 c 167 c 168 print*,' called sgeco - condition = ',1.0/rcond 169 c 170 c the test below might seem impossible to satisfy, but we'll see 171 c this is not so. this will be out jumping off point to discuss 172 c the availabe precision due to the finite word size for any computer 173 if (1.0 + rcond .le. 1.) then 174 print*,' Numerically singular matrix ' 175 c 176 c in a normal problem solving situation you would "stop" when a 177 c numerical singular matrix is detected. for cs 575, we will 178 c proceed to gain insight into the effects on accuracy 179 c 180 c a large condition number (effectively a small rcond = 1/cond) 181 c indicates extreme sensitivity in the matrix A - user beware! 182 c 183 endif 184 c ------------------------------------- 185 c using the LU factors of A (computed above by sgeco) and the 186 c RHS (stored in X initially) - compute the solution for A X = RHS 187 c 188 call sgesl (a, ndim, N, ipvt, X, 0) 189 c 190 print*,' the solution is x = ' 191 do row= k,1,-1 192 print *, ' row =', row 193 print 10, (X((row-1)*M + i), i=1,M) 194 end do 195 10 format (10f8.4) 196 stop f_diffusion.f Wed Oct 16 07:59:49 2002 page 4 197 end sgeco.f Wed Oct 16 07:59:49 2002 page 5 FILE "sgeco.f" 1 c *** from netlib, Wed Jan 30 12:37:38 EST 1991 *** 2 subroutine sgeco(a,lda,n,ipvt,rcond,z) ^ **** WAR #391: variable "lda" used only in unreached code 3 integer lda,n,ipvt(1) 4 real a(lda,1),z(1) 5 real rcond 6 c 7 c sgeco factors a real matrix by gaussian elimination 8 c and estimates the condition of the matrix. 9 c 10 c if rcond is not needed, sgefa is slightly faster. 11 c to solve a*x = b , follow sgeco by sgesl. 12 c to compute inverse(a)*c , follow sgeco by sgesl. 13 c to compute determinant(a) , follow sgeco by sgedi. 14 c to compute inverse(a) , follow sgeco by sgedi. 15 c 16 c on entry 17 c 18 c a real(lda, n) 19 c the matrix to be factored. 20 c 21 c lda integer 22 c the leading dimension of the array a . 23 c 24 c n integer 25 c the order of the matrix a . 26 c 27 c on return 28 c 29 c a an upper triangular matrix and the multipliers 30 c which were used to obtain it. 31 c the factorization can be written a = l*u where 32 c l is a product of permutation and unit lower 33 c triangular matrices and u is upper triangular. 34 c 35 c ipvt integer(n) 36 c an integer vector of pivot indices. 37 c 38 c rcond real 39 c an estimate of the reciprocal condition of a . 40 c for the system a*x = b , relative perturbations 41 c in a and b of size epsilon may cause 42 c relative perturbations in x of size epsilon/rcond . 43 c if rcond is so small that the logical expression 44 c 1.0 + rcond .eq. 1.0 45 c is true, then a may be singular to working 46 c precision. in particular, rcond is zero if 47 c exact singularity is detected or the estimate 48 c underflows. 49 c 50 c z real(n) 51 c a work vector whose contents are usually unimportant. 52 c if a is close to a singular matrix, then z is 53 c an approximate null vector in the sense that 54 c norm(a*z) = rcond*norm(a)*norm(z) . 55 c 56 c linpack. this version dated 08/14/78 . 57 c cleve moler, university of new mexico, argonne national lab. 58 c 59 c subroutines and functions 60 c 61 c linpack sgefa 62 c blas saxpy,sdot,sscal,sasum 63 c fortran abs,amax1,sign sgeco.f Wed Oct 16 07:59:49 2002 page 6 64 c 65 c internal variables 66 c 67 real sdot,ek,t,wk,wkm 68 real anorm,s,sasum,sm,ynorm 69 integer info,j,k,kb,kp1,l 70 c 71 c 72 c compute 1-norm of a 73 c 74 anorm = 0.0e0 75 do 10 j = 1, n 76 anorm = amax1(anorm,sasum(n,a(1,j),1)) 77 10 continue 78 c 79 c factor 80 c 81 call sgefa(a,lda,n,ipvt,info) ^ **** WAR #320: variable "info" set but never referenced 82 c 83 c rcond = 1/(norm(a)*(estimate of norm(inverse(a)))) . 84 c estimate = norm(z)/norm(y) where a*z = y and trans(a)*y = e . 85 c trans(a) is the transpose of a . the components of e are 86 c chosen to cause maximum local growth in the elements of w where 87 c trans(u)*w = e . the vectors are frequently rescaled to avoid 88 c overflow. 89 c 90 c solve trans(u)*w = e 91 c 92 ek = 1.0e0 93 do 20 j = 1, n 94 z(j) = 0.0e0 95 20 continue 96 do 100 k = 1, n 97 if (z(k) .ne. 0.0e0) ek = sign(ek,-z(k)) 98 if (abs(ek-z(k)) .le. abs(a(k,k))) go to 30 99 s = abs(a(k,k))/abs(ek-z(k)) 100 call sscal(n,s,z,1) 101 ek = s*ek 102 30 continue 103 wk = ek - z(k) 104 wkm = -ek - z(k) 105 s = abs(wk) 106 sm = abs(wkm) 107 if (a(k,k) .eq. 0.0e0) go to 40 108 wk = wk/a(k,k) 109 wkm = wkm/a(k,k) 110 go to 50 111 40 continue 112 wk = 1.0e0 113 wkm = 1.0e0 114 50 continue 115 kp1 = k + 1 116 if (kp1 .gt. n) go to 90 117 do 60 j = kp1, n 118 sm = sm + abs(z(j)+wkm*a(k,j)) 119 z(j) = z(j) + wk*a(k,j) 120 s = s + abs(z(j)) 121 60 continue 122 if (s .ge. sm) go to 80 123 t = wkm - wk 124 wk = wkm 125 do 70 j = kp1, n 126 z(j) = z(j) + t*a(k,j) 127 70 continue sgeco.f Wed Oct 16 07:59:49 2002 page 7 128 80 continue 129 90 continue 130 z(k) = wk 131 100 continue 132 s = 1.0e0/sasum(n,z,1) 133 call sscal(n,s,z,1) 134 c 135 c solve trans(l)*y = w 136 c 137 do 120 kb = 1, n 138 k = n + 1 - kb 139 if (k .lt. n) z(k) = z(k) + sdot(n-k,a(k+1,k),1,z(k+1),1) 140 if (abs(z(k)) .le. 1.0e0) go to 110 141 s = 1.0e0/abs(z(k)) 142 call sscal(n,s,z,1) 143 110 continue 144 l = ipvt(k) 145 t = z(l) 146 z(l) = z(k) 147 z(k) = t 148 120 continue 149 s = 1.0e0/sasum(n,z,1) 150 call sscal(n,s,z,1) 151 c 152 ynorm = 1.0e0 153 c 154 c solve l*v = y 155 c 156 do 140 k = 1, n 157 l = ipvt(k) 158 t = z(l) 159 z(l) = z(k) 160 z(k) = t 161 if (k .lt. n) call saxpy(n-k,t,a(k+1,k),1,z(k+1),1) 162 if (abs(z(k)) .le. 1.0e0) go to 130 163 s = 1.0e0/abs(z(k)) 164 call sscal(n,s,z,1) 165 ynorm = s*ynorm 166 130 continue 167 140 continue 168 s = 1.0e0/sasum(n,z,1) 169 call sscal(n,s,z,1) 170 ynorm = s*ynorm 171 c 172 c solve u*z = v 173 c 174 do 160 kb = 1, n 175 k = n + 1 - kb 176 if (abs(z(k)) .le. abs(a(k,k))) go to 150 177 s = abs(a(k,k))/abs(z(k)) 178 call sscal(n,s,z,1) 179 ynorm = s*ynorm 180 150 continue 181 if (a(k,k) .ne. 0.0e0) z(k) = z(k)/a(k,k) 182 if (a(k,k) .eq. 0.0e0) z(k) = 1.0e0 183 t = -z(k) 184 call saxpy(k-1,t,a(1,k),1,z(1),1) 185 160 continue 186 c make znorm = 1.0 187 s = 1.0e0/sasum(n,z,1) 188 call sscal(n,s,z,1) 189 ynorm = s*ynorm 190 c 191 if (anorm .ne. 0.0e0) rcond = ynorm/anorm 192 if (anorm .eq. 0.0e0) rcond = 0.0e0 193 return sgeco.f Wed Oct 16 07:59:49 2002 page 8 194 end 195 c============================================ 196 subroutine sgefa(a,lda,n,ipvt,info) ^ **** WAR #391: variable "lda" used only in unreached code 197 integer lda,n,ipvt(1),info 198 real a(lda,1) 199 c 200 c sgefa factors a real matrix by gaussian elimination. 201 c 202 c sgefa is usually called by sgeco, but it can be called 203 c directly with a saving in time if rcond is not needed. 204 c (time for sgeco) = (1 + 9/n)*(time for sgefa) . 205 c 206 c on entry 207 c 208 c a real(lda, n) 209 c the matrix to be factored. 210 c 211 c lda integer 212 c the leading dimension of the array a . 213 c 214 c n integer 215 c the order of the matrix a . 216 c 217 c on return 218 c 219 c a an upper triangular matrix and the multipliers 220 c which were used to obtain it. 221 c the factorization can be written a = l*u where 222 c l is a product of permutation and unit lower 223 c triangular matrices and u is upper triangular. 224 c 225 c ipvt integer(n) 226 c an integer vector of pivot indices. 227 c 228 c info integer 229 c = 0 normal value. 230 c = k if u(k,k) .eq. 0.0 . this is not an error 231 c condition for this subroutine, but it does 232 c indicate that sgesl or sgedi will divide by zero 233 c if called. use rcond in sgeco for a reliable 234 c indication of singularity. 235 c 236 c linpack. this version dated 08/14/78 . 237 c cleve moler, university of new mexico, argonne national lab. 238 c 239 c subroutines and functions 240 c 241 c blas saxpy,sscal,isamax 242 c 243 c internal variables 244 c 245 real t 246 integer isamax,j,k,kp1,l,nm1 247 c 248 c 249 c gaussian elimination with partial pivoting 250 c 251 info = 0 252 nm1 = n - 1 253 if (nm1 .lt. 1) go to 70 254 do 60 k = 1, nm1 255 kp1 = k + 1 256 c 257 c find l = pivot index sgeco.f Wed Oct 16 07:59:49 2002 page 9 258 c 259 l = isamax(n-k+1,a(k,k),1) + k - 1 260 ipvt(k) = l 261 c 262 c zero pivot implies this column already triangularized 263 c 264 if (a(l,k) .eq. 0.0e0) go to 40 265 c 266 c interchange if necessary 267 c 268 if (l .eq. k) go to 10 269 t = a(l,k) 270 a(l,k) = a(k,k) 271 a(k,k) = t 272 10 continue 273 c 274 c compute multipliers 275 c 276 t = -1.0e0/a(k,k) 277 call sscal(n-k,t,a(k+1,k),1) 278 c 279 c row elimination with column indexing 280 c 281 do 30 j = kp1, n 282 t = a(l,j) 283 if (l .eq. k) go to 20 284 a(l,j) = a(k,j) 285 a(k,j) = t 286 20 continue 287 call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 288 30 continue 289 go to 50 290 40 continue 291 info = k 292 50 continue 293 60 continue 294 70 continue 295 ipvt(n) = n 296 if (a(n,n) .eq. 0.0e0) info = n 297 return 298 end sgesl.f Wed Oct 16 07:59:49 2002 page 10 FILE "sgesl.f" 1 c *** from netlib, Mon Feb 4 14:58:46 EST 1991 *** 2 subroutine sgesl(a,lda,n,ipvt,b,job) ^ **** WAR #391: variable "lda" used only in unreached code 3 integer lda,n,ipvt(1),job 4 real a(lda,1),b(1) 5 c 6 c sgesl solves the real system 7 c a * x = b or trans(a) * x = b 8 c using the factors computed by sgeco or sgefa. 9 c 10 c on entry 11 c 12 c a real(lda, n) 13 c the output from sgeco or sgefa. 14 c 15 c lda integer 16 c the leading dimension of the array a . 17 c 18 c n integer 19 c the order of the matrix a . 20 c 21 c ipvt integer(n) 22 c the pivot vector from sgeco or sgefa. 23 c 24 c b real(n) 25 c the right hand side vector. 26 c 27 c job integer 28 c = 0 to solve a*x = b , 29 c = nonzero to solve trans(a)*x = b where 30 c trans(a) is the transpose. 31 c 32 c on return 33 c 34 c b the solution vector x . 35 c 36 c error condition 37 c 38 c a division by zero will occur if the input factor contains a 39 c zero on the diagonal. technically this indicates singularity 40 c but it is often caused by improper arguments or improper 41 c setting of lda . it will not occur if the subroutines are 42 c called correctly and if sgeco has set rcond .gt. 0.0 43 c or sgefa has set info .eq. 0 . 44 c 45 c to compute inverse(a) * c where c is a matrix 46 c with p columns 47 c call sgeco(a,lda,n,ipvt,rcond,z) 48 c if (rcond is too small) go to ... 49 c do 10 j = 1, p 50 c call sgesl(a,lda,n,ipvt,c(1,j),0) 51 c 10 continue 52 c 53 c linpack. this version dated 08/14/78 . 54 c cleve moler, university of new mexico, argonne national lab. 55 c 56 c subroutines and functions 57 c 58 c blas saxpy,sdot 59 c 60 c internal variables 61 c 62 real sdot,t 63 integer k,kb,l,nm1 sgesl.f Wed Oct 16 07:59:49 2002 page 11 64 c 65 nm1 = n - 1 66 if (job .ne. 0) go to 50 67 c 68 c job = 0 , solve a * x = b 69 c first solve l*y = b 70 c 71 if (nm1 .lt. 1) go to 30 72 do 20 k = 1, nm1 73 l = ipvt(k) 74 t = b(l) 75 if (l .eq. k) go to 10 76 b(l) = b(k) 77 b(k) = t 78 10 continue 79 call saxpy(n-k,t,a(k+1,k),1,b(k+1),1) 80 20 continue 81 30 continue 82 c 83 c now solve u*x = y 84 c 85 do 40 kb = 1, n 86 k = n + 1 - kb 87 b(k) = b(k)/a(k,k) 88 t = -b(k) 89 call saxpy(k-1,t,a(1,k),1,b(1),1) 90 40 continue 91 go to 100 92 50 continue 93 c 94 c job = nonzero, solve trans(a) * x = b 95 c first solve trans(u)*y = b 96 c 97 do 60 k = 1, n 98 t = sdot(k-1,a(1,k),1,b(1),1) 99 b(k) = (b(k) - t)/a(k,k) 100 60 continue 101 c 102 c now solve trans(l)*x = y 103 c 104 if (nm1 .lt. 1) go to 90 105 do 80 kb = 1, nm1 106 k = n - kb 107 b(k) = b(k) + sdot(n-k,a(k+1,k),1,b(k+1),1) 108 l = ipvt(k) 109 if (l .eq. k) go to 70 110 t = b(l) 111 b(l) = b(k) 112 b(k) = t 113 70 continue 114 80 continue 115 90 continue 116 100 continue 117 return 118 end 119 blas.f Wed Oct 16 07:59:49 2002 page 12 FILE "blas.f" 1 subroutine sscal(n,sa,sx,incx) 2 c 3 c scales a vector by a constant. 4 c uses unrolled loops for increment equal to 1. 5 c jack dongarra, linpack, 3/11/78. 6 c 7 real sa,sx(1) 8 integer i,incx,m,mp1,n,nincx 9 c 10 if(n.le.0)return 11 if(incx.eq.1)go to 20 12 c 13 c code for increment not equal to 1 14 c 15 nincx = n*incx 16 do 10 i = 1,nincx,incx 17 sx(i) = sa*sx(i) 18 10 continue 19 return 20 c 21 c code for increment equal to 1 22 c 23 c 24 c clean-up loop 25 c 26 20 m = mod(n,5) 27 if( m .eq. 0 ) go to 40 28 do 30 i = 1,m 29 sx(i) = sa*sx(i) 30 30 continue 31 if( n .lt. 5 ) return 32 40 mp1 = m + 1 33 do 50 i = mp1,n,5 34 sx(i) = sa*sx(i) 35 sx(i + 1) = sa*sx(i + 1) 36 sx(i + 2) = sa*sx(i + 2) 37 sx(i + 3) = sa*sx(i + 3) 38 sx(i + 4) = sa*sx(i + 4) 39 50 continue 40 return 41 end 42 integer function isamax(n,sx,incx) 43 c 44 c finds the index of element having max. absolute value. 45 c jack dongarra, linpack, 3/11/78. 46 c 47 real sx(1),smax 48 integer i,incx,ix,n 49 c 50 isamax = 0 51 if( n .lt. 1 ) return 52 isamax = 1 53 if(n.eq.1)return 54 if(incx.eq.1)go to 20 55 c 56 c code for increment not equal to 1 57 c 58 ix = 1 59 smax = abs(sx(1)) 60 ix = ix + incx 61 do 10 i = 2,n 62 if(abs(sx(ix)).le.smax) go to 5 63 isamax = i 64 smax = abs(sx(ix)) 65 5 ix = ix + incx blas.f Wed Oct 16 07:59:49 2002 page 13 66 10 continue 67 return 68 c 69 c code for increment equal to 1 70 c 71 20 smax = abs(sx(1)) 72 do 30 i = 2,n 73 if(abs(sx(i)).le.smax) go to 30 74 isamax = i 75 smax = abs(sx(i)) 76 30 continue 77 return 78 end 79 real function sasum(n,sx,incx) 80 c 81 c takes the sum of the absolute values. 82 c uses unrolled loops for increment equal to one. 83 c jack dongarra, linpack, 3/11/78. 84 c 85 real sx(1),stemp 86 integer i,incx,m,mp1,n,nincx 87 c 88 sasum = 0.0e0 89 stemp = 0.0e0 90 if(n.le.0)return 91 if(incx.eq.1)go to 20 92 c 93 c code for increment not equal to 1 94 c 95 nincx = n*incx 96 do 10 i = 1,nincx,incx 97 stemp = stemp + abs(sx(i)) 98 10 continue 99 sasum = stemp 100 return 101 c 102 c code for increment equal to 1 103 c 104 c 105 c clean-up loop 106 c 107 20 m = mod(n,6) 108 if( m .eq. 0 ) go to 40 109 do 30 i = 1,m 110 stemp = stemp + abs(sx(i)) 111 30 continue 112 if( n .lt. 6 ) go to 60 113 40 mp1 = m + 1 114 do 50 i = mp1,n,6 115 stemp = stemp + abs(sx(i)) + abs(sx(i + 1)) + abs(sx(i + 2)) 116 * + abs(sx(i + 3)) + abs(sx(i + 4)) + abs(sx(i + 5)) 117 50 continue 118 60 sasum = stemp 119 return 120 end 121 subroutine saxpy(n,sa,sx,incx,sy,incy) 122 c 123 c constant times a vector plus a vector. 124 c uses unrolled loop for increments equal to one. 125 c jack dongarra, linpack, 3/11/78. 126 c 127 real sx(1),sy(1),sa 128 integer i,incx,incy,ix,iy,m,mp1,n 129 c 130 if(n.le.0)return 131 if (sa .eq. 0.0) return blas.f Wed Oct 16 07:59:49 2002 page 14 132 if(incx.eq.1.and.incy.eq.1)go to 20 133 c 134 c code for unequal increments or equal increments 135 c not equal to 1 136 c 137 ix = 1 138 iy = 1 139 if(incx.lt.0)ix = (-n+1)*incx + 1 140 if(incy.lt.0)iy = (-n+1)*incy + 1 141 do 10 i = 1,n 142 sy(iy) = sy(iy) + sa*sx(ix) 143 ix = ix + incx 144 iy = iy + incy 145 10 continue 146 return 147 c 148 c code for both increments equal to 1 149 c 150 c 151 c clean-up loop 152 c 153 20 m = mod(n,4) 154 if( m .eq. 0 ) go to 40 155 do 30 i = 1,m 156 sy(i) = sy(i) + sa*sx(i) 157 30 continue 158 if( n .lt. 4 ) return 159 40 mp1 = m + 1 160 do 50 i = mp1,n,4 161 sy(i) = sy(i) + sa*sx(i) 162 sy(i + 1) = sy(i + 1) + sa*sx(i + 1) 163 sy(i + 2) = sy(i + 2) + sa*sx(i + 2) 164 sy(i + 3) = sy(i + 3) + sa*sx(i + 3) 165 50 continue 166 return 167 end 168 real function sdot(n,sx,incx,sy,incy) 169 c 170 c forms the dot product of two vectors. 171 c uses unrolled loops for increments equal to one. 172 c jack dongarra, linpack, 3/11/78. 173 c 174 real sx(1),sy(1),stemp 175 integer i,incx,incy,ix,iy,m,mp1,n 176 c 177 stemp = 0.0e0 178 sdot = 0.0e0 179 if(n.le.0)return 180 if(incx.eq.1.and.incy.eq.1)go to 20 181 c 182 c code for unequal increments or equal increments 183 c not equal to 1 184 c 185 ix = 1 186 iy = 1 187 if(incx.lt.0)ix = (-n+1)*incx + 1 188 if(incy.lt.0)iy = (-n+1)*incy + 1 189 do 10 i = 1,n 190 stemp = stemp + sx(ix)*sy(iy) 191 ix = ix + incx 192 iy = iy + incy 193 10 continue 194 sdot = stemp 195 return 196 c 197 c code for both increments equal to 1 blas.f Wed Oct 16 07:59:49 2002 page 15 198 c 199 c 200 c clean-up loop 201 c 202 20 m = mod(n,5) 203 if( m .eq. 0 ) go to 40 204 do 30 i = 1,m 205 stemp = stemp + sx(i)*sy(i) 206 30 continue 207 if( n .lt. 5 ) go to 60 208 40 mp1 = m + 1 209 do 50 i = mp1,n,5 210 stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + 211 * sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4 ) 212 50 continue 213 60 sdot = stemp 214 return 215 end 216 217 Cross Reference Wed Oct 16 07:59:49 2002 page 16 C R O S S R E F E R E N C E T A B L E Source files: f_diffusion.f sgeco.f sgesl.f blas.f Legend: D Definition/Declaration U Simple use M Modified occurrence A Actual argument C Subroutine/Function call I Initialization: DATA or extended declaration E Occurrence in EQUIVALENCE N Occurrence in NAMELIST L Use Module Cross Reference Wed Oct 16 07:59:49 2002 page 17 P R O G R A M F O R M Program ------- f_diffusion f_diffusion.f D 1:D Cross Reference Wed Oct 16 07:59:49 2002 page 18 Functions and Subroutines ------------------------- ABS intrinsic sgeco.f C 98:C 98:C 99:C 99:C 105:C 106:C 118:C 120:C 140:C 141:C 162:C 163:C 176:C 176:C 177:C 177:C blas.f C 59:C 62:C 64:C 71:C 73:C 75:C C 97:C 110:C 115:C 115:C 115:C 116:C 116:C 116:C AMAX1 intrinsic sgeco.f C 76:C SIGN intrinsic sgeco.f C 97:C isamax int*4 sgeco.f DC 246:D 259:C blas.f DM 42:D 50:M 52:M 63:M 74:M max_r_rr intrinsic sgeco.f C 4:C C 198:C sgesl.f C 4:C sasum real*4 sgeco.f DC 68:D 76:C 132:C 149:C 168:C 187:C blas.f DM 79:D 88:M 99:M 118:M saxpy sgeco.f C 161:C 184:C C 287:C sgesl.f C 79:C 89:C blas.f D 121:D sdot real*4 sgeco.f DC 67:D 139:C sgesl.f DC 62:D 98:C 107:C blas.f DM 168:D 178:M 194:M 213:M sgeco f_diffusion.f C 165:C sgeco.f D 2:D sgefa sgeco.f C 81:C D 196:D sgesl f_diffusion.f C 188:C sgesl.f D 2:D sscal sgeco.f C 100:C 133:C 142:C 150:C 164:C 169:C 178:C 188:C C 277:C blas.f D 1:D Cross Reference Wed Oct 16 07:59:49 2002 page 19 Parameters ---------- ndim int*4 f_diffusion.f DUA 2:D 44:U 46:U 165:A 188:A sidebottom real*4 f_diffusion.f DU 3:D 127:U 130:U sideleft real*4 f_diffusion.f DU 2:D 127:U 138:U 145:U sideright real*4 f_diffusion.f DU 2:D 133:U 141:U 153:U sidetop real*4 f_diffusion.f DU 3:D 133:U 145:U 149:U 153:U Cross Reference Wed Oct 16 07:59:49 2002 page 20 Variables and Arrays -------------------- a real*4 array size: 16777216 bytes f_diffusion.f DMA 22:D 54:M 64:M 65:M 66:M 69:M 70:M 71:M 72:M 75:M 76:M 77:M 86:M 87:M 88:M 89:M 92:M 93:M 94:M 95:M 96:M 100:M 101:M 102:M 103:M 108:M 109:M 110:M 113:M 114:M 115:M 116:M 119:M 120:M 121:M 165:A 188:A a real*4 dummy array sgeco.f DUA 2:D 4:D 76:A 81:A 98:A 99:A 107:U 108:U 109:U 118:A 119:U 126:U 139:A 161:A 173:U 176:A 177:A 181:U 181:U 182:U 184:A sgesl.f DUA 2:D 4:D 79:A 87:U 89:A 98:A 99:U 104:U 107:A sgeco.f DUMA 196:D 198:D 259:A 264:U 269:U 270:M 270:U 271:M 276:U 277:A 280:U 282:U 284:M 284:U 285:M 287:A 287:A 296:U anorm real*4 sgeco.f DUMA 68:D 74:M 76:M 76:A 191:U 191:U 192:U b real*4 dummy array size: 4 bytes sgesl.f DUMA 2:D 4:D 74:U 76:M 76:U 77:M 79:A 87:M 87:U 88:U 89:A 98:A 99:M 99:U 107:M 107:U 107:A 110:U 111:M 111:U 112:M col int*4 f_diffusion.f DM 31:D 90:M 111:M 146:M ek real*4 sgeco.f DUMA 67:D 92:M 97:M 97:A 98:A 99:A 101:M 101:U 103:U 104:U 104:U i int*4 f_diffusion.f UM 67:M 193:U blas.f DUM 86:D 96:M 97:U 109:M 110:U 114:M 115:U 115:U 115:U 116:U 116:U 116:U DUM 175:D 189:M 204:M 205:U 205:U 209:M 210:U 210:U 210:U 210:U 211:U 211:U 211:U 211:U 211:U 211:U DUM 8:D 16:M 17:U 17:U 28:M 29:U 29:U 33:M 34:U 34:U 35:U 35:U 36:U 36:U 37:U 37:U 38:U 38:U DUM 128:D 141:M 155:M 156:U 156:U 156:U 160:M 161:U 161:U 161:U 162:U 162:U 162:U 163:U 163:U 163:U 164:U 164:U 164:U DUM 48:D 61:M 63:U 72:M Cross Reference Wed Oct 16 07:59:49 2002 page 21 73:U 74:U 75:U incx int*4 dummy blas.f DU 79:D 86:D 91:U 95:U 96:U DU 168:D 175:D 180:U 187:U 187:U 191:U DU 1:D 8:D 11:U 15:U 16:U DU 121:D 128:D 132:U 139:U 139:U 143:U DU 42:D 48:D 54:U 60:U 65:U incy int*4 dummy blas.f DU 168:D 175:D 180:U 188:U 188:U 192:U DU 121:D 128:D 132:U 140:U 140:U 144:U index int*4 f_diffusion.f UM 63:M 64:U 64:U 65:U 65:U 66:U 66:U 68:M 68:U 69:U 69:U 70:U 70:U 71:U 71:U 72:U 72:U 74:M 74:U 75:U 75:U 76:U 76:U 77:U 77:U 84:M 84:U 86:U 86:U 87:U 87:U 88:U 88:U 89:U 89:U 91:M 91:U 92:U 92:U 93:U 93:U 94:U 94:U 95:U 95:U 96:U 96:U 99:M 99:U 100:U 100:U 101:U 101:U 102:U 102:U 103:U 103:U 107:M 107:U 108:U 108:U 109:U 109:U 110:U 110:U 112:M 112:U 113:U 113:U 114:U 114:U 115:U 115:U 116:U 116:U 118:M 118:U 119:U 119:U 120:U 120:U 121:U 121:U 128:M 130:U 134:M 136:M 136:U 138:U 139:M 139:U 141:U 144:M 144:U 145:U 147:M 147:U 149:U 152:M 152:U 153:U info int*4 sgeco.f DA 69:D 81:A info int*4 dummy sgeco.f DM 196:D 197:D 251:M 291:M 296:M ipvt int*4 array size: 8192 bytes f_diffusion.f DA 31:D 165:A 188:A ipvt int*4 dummy array size: 4 bytes Cross Reference Wed Oct 16 07:59:49 2002 page 22 sgeco.f DUA 2:D 3:D 81:A 144:U 157:U sgesl.f DU 2:D 3:D 73:U 108:U sgeco.f DM 196:D 197:D 260:M 295:M ix int*4 blas.f DUM 175:D 185:M 187:M 190:U 191:M 191:U DUM 128:D 137:M 139:M 142:U 143:M 143:U DUM 48:D 58:M 60:M 60:U 62:U 64:U 65:M 65:U iy int*4 blas.f DUM 175:D 186:M 188:M 190:U 192:M 192:U DUM 128:D 138:M 140:M 142:U 142:U 144:M 144:U j int*4 sgeco.f DUM 69:D 75:M 76:U 93:M 94:U 117:M 118:U 118:U 119:U 119:U 119:U 120:U 125:M 126:U 126:U 126:U DUM 246:D 281:M 282:U 284:U 284:U 285:U 287:U job int*4 dummy sgesl.f DU 2:D 3:D 66:U k int*4 f_diffusion.f UM 39:M 43:U 83:U 135:U 191:U sgeco.f DUM 69:D 96:M 97:U 97:U 98:U 98:U 98:U 99:U 99:U 99:U 103:U 104:U 107:U 107:U 108:U 108:U 109:U 109:U 115:U 118:U 119:U 126:U 130:U 138:M 139:U 139:U 139:U 139:U 139:U 139:U 139:U 140:U 141:U 144:U 146:U 147:U 156:M 157:U 159:U 160:U 161:U 161:U 161:U 161:U 161:U 162:U 163:U 175:M 176:U 176:U 176:U 177:U 177:U 177:U 181:U 181:U 181:U 181:U 181:U 181:U 182:U 182:U 182:U 183:U 184:U 184:U sgesl.f DUM 63:D 72:M 73:U 75:U 76:U 77:U 79:U 79:U 79:U 79:U 86:M 87:U 87:U 87:U 87:U 88:U 89:U 89:U 97:M 98:U 98:U 99:U 99:U 99:U 99:U 106:M 107:U 107:U 107:U 107:U 107:U 107:U 108:U 109:U 111:U 112:U sgeco.f DUM 246:D 254:M 255:U 259:U 259:U 259:U 259:U 260:U 264:U 268:U 269:U 270:U 270:U 270:U 271:U 271:U 276:U 276:U 277:U 277:U 277:U 283:U 284:U 285:U 287:U 287:U 287:U 287:U 291:U kb int*4 sgeco.f DUM 69:D 137:M 138:U 174:M 175:U sgesl.f DUM 63:D 85:M 86:U 105:M 106:U kp1 int*4 sgeco.f DUM 69:D 115:M 116:U 117:U 125:U DUM 246:D 255:M 281:U Cross Reference Wed Oct 16 07:59:49 2002 page 23 l int*4 sgeco.f DUM 69:D 144:M 145:U 146:U 157:M 158:U 159:U sgesl.f DUM 63:D 73:M 74:U 75:U 76:U 108:M 109:U 110:U 111:U sgeco.f DUM 246:D 259:M 260:U 264:U 268:U 269:U 270:U 282:U 283:U 284:U lda int*4 dummy sgeco.f DUA 2:D 3:D 4:U 81:A sgesl.f DU 2:D 3:D 4:U sgeco.f DU 196:D 197:D 198:U m int*4 f_diffusion.f UM 42:M 43:U 66:U 67:U 72:U 77:U 88:U 89:U 90:U 95:U 96:U 102:U 103:U 109:U 111:U 116:U 121:U 128:U 133:U 134:U 139:U 146:U 193:U blas.f DUM 86:D 107:M 108:U 109:U 113:U DUM 175:D 202:M 203:U 204:U 208:U DUM 8:D 26:M 27:U 28:U 32:U DUM 128:D 153:M 154:U 155:U 159:U mp1 int*4 blas.f DUM 86:D 113:M 114:U DUM 175:D 208:M 209:U DUM 8:D 32:M 33:U DUM 128:D 159:M 160:U n int*4 f_diffusion.f UMA 43:M 44:U 54:U 54:U 125:U 165:A 188:A n int*4 dummy sgeco.f DUA 2:D 3:D 75:U 76:A 81:A 93:U 96:U 100:A 116:U 117:U 125:U 132:A 133:A 137:U 138:U 139:U 139:U 142:A 149:A 150:A 156:U 161:U 161:U 164:A 168:A 169:A 174:U 175:U 178:A 187:A 188:A sgesl.f DU 2:D 3:D 65:U 79:U 85:U 86:U 97:U 106:U 107:U blas.f DU 79:D 86:D 90:U 95:U 107:U 112:U 114:U DU 168:D 175:D 179:U 187:U 188:U 189:U 202:U 207:U 209:U sgeco.f DU 196:D 197:D 252:U 259:U 277:U 281:U 287:U 295:U 295:U 296:U 296:U 296:U blas.f DU 1:D 8:D 10:U 15:U 26:U 31:U 33:U DU 121:D 128:D 130:U 139:U 140:U 141:U 153:U 158:U 160:U DU 42:D 48:D 51:U 53:U 61:U 72:U nincx int*4 blas.f DUM 86:D 95:M 96:U DUM 8:D 15:M 16:U Cross Reference Wed Oct 16 07:59:49 2002 page 24 nm1 int*4 sgesl.f DUM 63:D 65:M 71:U 72:U 104:U 105:U sgeco.f DUM 246:D 252:M 253:U 254:U rcond real*4 f_diffusion.f UA 165:A 168:U 173:U rcond real*4 dummy sgeco.f DM 2:D 5:D 191:M 192:M row int*4 f_diffusion.f DUM 31:D 83:M 135:M 191:M 192:U 193:U s real*4 sgeco.f DUMA 68:D 99:M 100:A 101:U 105:M 120:M 120:U 122:U 132:M 133:A 141:M 142:A 149:M 150:A 163:M 164:A 165:U 168:M 169:A 170:U 177:M 178:A 179:U 187:M 188:A 189:U sa real*4 dummy blas.f DU 1:D 7:D 17:U 29:U 34:U 35:U 36:U 37:U 38:U DU 121:D 127:D 131:U 142:U 156:U 161:U 162:U 163:U 164:U sm real*4 sgeco.f DUM 68:D 106:M 118:M 118:U 122:U smax real*4 blas.f DUM 47:D 59:M 62:U 64:M 71:M 73:U 75:M stemp real*4 blas.f DUM 85:D 89:M 97:M 97:U 99:U 110:M 110:U 115:M 115:U 118:U DUM 174:D 177:M 190:M 190:U 194:U 205:M 205:U 210:M 210:U 213:U sx real*4 dummy array size: 4 bytes blas.f DA 79:D 85:D 97:A 110:A 115:A 115:A 115:A 116:A 116:A 116:A DU 168:D 174:D 190:U 205:U 210:U 210:U 211:U 211:U 211:U DUM 1:D 7:D 17:M 17:U 29:M 29:U 34:M 34:U 35:M 35:U 36:M 36:U 37:M 37:U 38:M 38:U DU 121:D 127:D 142:U 156:U 161:U 162:U 163:U 164:U DA 42:D 47:D 59:A 62:A 64:A 71:A 73:A 75:A sy real*4 dummy array size: 4 bytes blas.f DU 168:D 174:D 190:U 205:U 210:U 210:U 211:U 211:U 211:U DUM 121:D 127:D 142:M 142:U 156:M 156:U 161:M 161:U 162:M 162:U 163:M 163:U 164:M 164:U t real*4 sgeco.f DUMA 67:D 123:M 126:U 145:M 147:U 158:M 160:U 161:A 183:M 184:A Cross Reference Wed Oct 16 07:59:49 2002 page 25 sgesl.f DUMA 62:D 74:M 77:U 79:A 88:M 89:A 98:M 99:U 110:M 112:U sgeco.f DUMA 245:D 269:M 271:U 276:M 277:A 282:M 285:U 287:A wk real*4 sgeco.f DUMA 67:D 103:M 105:A 108:M 108:U 112:M 119:U 123:U 124:M 130:U wkm real*4 sgeco.f DUMA 67:D 104:M 106:A 109:M 109:U 113:M 118:A 123:U 124:U work real*4 array size: 8192 bytes f_diffusion.f DA 28:D 165:A x real*4 array size: 8192 bytes f_diffusion.f DUMA 25:D 125:M 127:M 130:M 133:M 138:M 141:M 145:M 149:M 153:M 188:A 193:U ynorm real*4 sgeco.f DUM 68:D 152:M 165:M 165:U 170:M 170:U 179:M 179:U 189:M 189:U 191:U z real*4 dummy array size: 4 bytes sgeco.f DUMA 2:D 4:D 94:M 97:U 97:A 98:A 99:A 100:A 103:U 104:U 118:A 119:M 119:U 120:A 126:M 126:U 130:M 132:A 133:A 139:M 139:U 139:A 140:A 141:A 142:A 145:U 146:M 146:U 147:M 149:A 150:A 158:U 159:M 159:U 160:M 161:A 162:A 163:A 164:A 168:A 169:A 176:A 177:A 178:A 181:M 181:U 182:M 183:U 184:A 187:A 188:A Cross Reference Wed Oct 16 07:59:49 2002 page 26 Format and Executable Labels ---------------------------- 5 blas.f DU 62:U 65:D 10 format f_diffusion.f DU 193:U 195:D 10 sgeco.f DU 75:U 77:D sgesl.f DU 75:U 78:D blas.f DU 96:U 98:D DU 189:U 193:D sgeco.f DU 268:U 272:D blas.f DU 16:U 18:D DU 141:U 145:D DU 61:U 66:D 20 sgeco.f DU 93:U 95:D sgesl.f DU 72:U 80:D blas.f DU 91:U 107:D DU 180:U 202:D sgeco.f DU 283:U 286:D blas.f DU 11:U 26:D DU 132:U 153:D DU 54:U 71:D 30 sgeco.f DU 98:U 102:D sgesl.f DU 71:U 81:D blas.f DU 109:U 111:D DU 204:U 206:D sgeco.f DU 281:U 288:D blas.f DU 28:U 30:D DU 155:U 157:D DU 73:U 76:D 40 sgeco.f DU 107:U 111:D sgesl.f DU 85:U 90:D blas.f DU 108:U 113:D DU 203:U 208:D sgeco.f DU 264:U 290:D blas.f DU 27:U 32:D DU 154:U 159:D 50 sgeco.f DU 110:U 114:D sgesl.f DU 66:U 92:D blas.f DU 114:U 117:D DU 209:U 212:D sgeco.f DU 227:U 292:D blas.f DU 33:U 39:D DU 160:U 165:D 60 sgeco.f DU 117:U 121:D sgesl.f DU 97:U 100:D blas.f DU 112:U 118:D DU 207:U 213:D sgeco.f DU 254:U 293:D 70 sgeco.f DU 125:U 127:D sgesl.f DU 109:U 113:D sgeco.f DU 253:U 294:D 80 sgeco.f DU 122:U 128:D sgesl.f DU 105:U 114:D 90 sgeco.f DU 116:U 129:D sgesl.f DU 104:U 115:D Cross Reference Wed Oct 16 07:59:49 2002 page 27 100 sgeco.f DU 96:U 131:D sgesl.f DU 91:U 116:D 110 sgeco.f DU 140:U 143:D 120 sgeco.f DU 137:U 148:D 130 sgeco.f DU 162:U 166:D 140 sgeco.f DU 156:U 167:D 150 sgeco.f DU 176:U 180:D 160 sgeco.f DU 174:U 185:D 99999 f_diffusion.f DU 36:D 47:U ------------------------------------------------------------------------------ STATISTIC Wed Oct 16 07:59:49 2002 page 28 Date: Wed Oct 16 07:59:49 2002 Options: /opt/SUNWspro/bin/../WS6U2/bin/f90 -Xlist Files: 4 Lines: 828 Routines: 9 (MAIN: 1; Subroutines: 5; Functions: 3) Messages: 4 (Warnings: 4)