| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | subroutine out_acsr(acsr,ia,ja,n)
|
|---|
| 3 | real acsr(1)
|
|---|
| 4 | integer ia(2),ja(1)
|
|---|
| 5 | character*1 adum
|
|---|
| 6 | c
|
|---|
| 7 | open (unit=40,file='q0')
|
|---|
| 8 | open (unit=41,file='q1')
|
|---|
| 9 | open (unit=42,file='q2')
|
|---|
| 10 | c
|
|---|
| 11 | write(6,*) 'this is ia:',ia(1),ia(2),n
|
|---|
| 12 | do i=1,n
|
|---|
| 13 | nc = ia(i+1)-ia(i)
|
|---|
| 14 | j0 = ia(i)-1
|
|---|
| 15 | do jc=1,nc
|
|---|
| 16 | j = ja (jc+j0)
|
|---|
| 17 | a = acsr(jc+j0)
|
|---|
| 18 | write(40,40) a
|
|---|
| 19 | write(41,41) i
|
|---|
| 20 | write(42,41) j
|
|---|
| 21 | write(6,*) i,j,a
|
|---|
| 22 | enddo
|
|---|
| 23 | enddo
|
|---|
| 24 | 40 format(1pe20.12)
|
|---|
| 25 | 41 format(i9)
|
|---|
| 26 | c
|
|---|
| 27 | close (unit=40)
|
|---|
| 28 | close (unit=41)
|
|---|
| 29 | close (unit=42)
|
|---|
| 30 | c
|
|---|
| 31 | c write(6,*) 'Output in q0,q1,q2. Continue? <cr> '
|
|---|
| 32 | c read (5,*,end=1,err=1) adum
|
|---|
| 33 | 1 continue
|
|---|
| 34 | c
|
|---|
| 35 | return
|
|---|
| 36 | end
|
|---|
| 37 | c-----------------------------------------------------------------------
|
|---|
| 38 | subroutine compress_acsr(acsr,ia,ja,n)
|
|---|
| 39 | c
|
|---|
| 40 | c Compress csr formatted a matrix based on drop tolerance, eps.
|
|---|
| 41 | c
|
|---|
| 42 | real acsr(1)
|
|---|
| 43 | integer ia(1),ja(1)
|
|---|
| 44 | c
|
|---|
| 45 | eps = 1.e-10
|
|---|
| 46 | c
|
|---|
| 47 | k0 = ia(1)-1
|
|---|
| 48 | do i=1,n
|
|---|
| 49 | nc = ia(i+1)-ia(i)
|
|---|
| 50 | j0 = ia(i)-1
|
|---|
| 51 |
|
|---|
| 52 | C row sum as metric
|
|---|
| 53 | aa = 0.0
|
|---|
| 54 | do jc=1,nc
|
|---|
| 55 | aa = aa+abs(acsr(jc+j0))
|
|---|
| 56 | enddo
|
|---|
| 57 | if (aa.lt.eps) then
|
|---|
| 58 | write(6,*) 'inf-norm of row is small',aa
|
|---|
| 59 | aa=1.0
|
|---|
| 60 | endif
|
|---|
| 61 | C
|
|---|
| 62 | kc = 0
|
|---|
| 63 | do jc=1,nc
|
|---|
| 64 | ae = abs(acsr(jc+j0)/aa)
|
|---|
| 65 | if (ae.gt.eps) then
|
|---|
| 66 | c bump column pointer
|
|---|
| 67 | kc = kc+1
|
|---|
| 68 | acsr(kc+k0) = acsr(jc+j0)
|
|---|
| 69 | ja (kc+k0) = ja (jc+j0)
|
|---|
| 70 | endif
|
|---|
| 71 | enddo
|
|---|
| 72 | ia (i) = k0+1
|
|---|
| 73 | k0 = k0+kc
|
|---|
| 74 | enddo
|
|---|
| 75 | ia(i) = k0+1
|
|---|
| 76 | c
|
|---|
| 77 | return
|
|---|
| 78 | end
|
|---|
| 79 | c-----------------------------------------------------------------------
|
|---|
| 80 | subroutine outbox(xmax,xmin,ymax,ymin,io)
|
|---|
| 81 | c
|
|---|
| 82 | dx = xmax-xmin
|
|---|
| 83 | dy = ymax-ymin
|
|---|
| 84 | dd = max(dx,dy)/2.
|
|---|
| 85 | c
|
|---|
| 86 | xh = (xmax+xmin)/2.
|
|---|
| 87 | yh = (ymax+ymin)/2.
|
|---|
| 88 | xmn = xh - dd
|
|---|
| 89 | xmx = xh + dd
|
|---|
| 90 | ymn = yh - dd
|
|---|
| 91 | ymx = yh + dd
|
|---|
| 92 | c
|
|---|
| 93 | write(io,*)
|
|---|
| 94 | write(io,*) ymn,xmn
|
|---|
| 95 | write(io,*) ymx,xmn
|
|---|
| 96 | write(io,*) ymx,xmx
|
|---|
| 97 | write(io,*) ymn,xmx
|
|---|
| 98 | write(io,*) ymn,xmn
|
|---|
| 99 | write(io,*)
|
|---|
| 100 | c
|
|---|
| 101 | return
|
|---|
| 102 | end
|
|---|
| 103 | c-----------------------------------------------------------------------
|
|---|
| 104 | subroutine imout(x,m,n,name)
|
|---|
| 105 | integer x(m,1)
|
|---|
| 106 | character*9 name
|
|---|
| 107 | c
|
|---|
| 108 | n15 = min(n,30)
|
|---|
| 109 | do i=1,m
|
|---|
| 110 | write(6,6) n,name,(x(i,k),k=1,n15)
|
|---|
| 111 | enddo
|
|---|
| 112 | 6 format(i4,1x,a9,':',2x,30i3)
|
|---|
| 113 | return
|
|---|
| 114 | end
|
|---|
| 115 | c-----------------------------------------------------------------------
|
|---|
| 116 | subroutine out_abd(abd,lda,n,m)
|
|---|
| 117 | real abd(lda,1)
|
|---|
| 118 | c
|
|---|
| 119 | write(6,*) 'abd:',lda,m,n
|
|---|
| 120 | do j=1,n
|
|---|
| 121 | write(6,6) (abd(i,j),i=lda,1,-1)
|
|---|
| 122 | enddo
|
|---|
| 123 | 6 format(12f8.3)
|
|---|
| 124 | c.
|
|---|
| 125 | c.
|
|---|
| 126 | open(unit=40,file='b0')
|
|---|
| 127 | open(unit=41,file='b1')
|
|---|
| 128 | open(unit=42,file='b2')
|
|---|
| 129 | mm = lda-1
|
|---|
| 130 | do 20 j = 1, n
|
|---|
| 131 | i1 = max0(1, j-mm)
|
|---|
| 132 | do 10 i = i1, j
|
|---|
| 133 | k = i-j+mm+1
|
|---|
| 134 | a = abd(k,j)
|
|---|
| 135 | write(40,*) a
|
|---|
| 136 | write(41,*) i
|
|---|
| 137 | write(42,*) j
|
|---|
| 138 | 10 continue
|
|---|
| 139 | 20 continue
|
|---|
| 140 | c
|
|---|
| 141 | close(unit=40)
|
|---|
| 142 | close(unit=41)
|
|---|
| 143 | close(unit=42)
|
|---|
| 144 | c
|
|---|
| 145 | return
|
|---|
| 146 | end
|
|---|
| 147 | c-----------------------------------------------------------------------
|
|---|
| 148 | subroutine rarr_out(x,name13)
|
|---|
| 149 | include 'SIZE'
|
|---|
| 150 | include 'INPUT'
|
|---|
| 151 | c
|
|---|
| 152 | real x(lx1,ly1,lz1,lelt)
|
|---|
| 153 | character*13 name13
|
|---|
| 154 | c
|
|---|
| 155 | if (nelv.gt.20) return
|
|---|
| 156 | write(6,*)
|
|---|
| 157 | write(6,1) name13
|
|---|
| 158 | 1 format('rarr ',3x,a13)
|
|---|
| 159 | c
|
|---|
| 160 | c 3 D
|
|---|
| 161 | c
|
|---|
| 162 | if (if3d) then
|
|---|
| 163 | do iz=1,lz1
|
|---|
| 164 | write(6,*)
|
|---|
| 165 | do j=ly1,1,-1
|
|---|
| 166 | write(6,3) (x(k,j,iz,1),k=1,lx1),(x(k,j,iz,2),k=1,lx1)
|
|---|
| 167 | enddo
|
|---|
| 168 | enddo
|
|---|
| 169 | write(6,*)
|
|---|
| 170 | write(6,*)
|
|---|
| 171 | return
|
|---|
| 172 | endif
|
|---|
| 173 | c
|
|---|
| 174 | c 2 D
|
|---|
| 175 | c
|
|---|
| 176 | if (nelv.gt.2) then
|
|---|
| 177 | write(6,*)
|
|---|
| 178 | do j=ly1,1,-1
|
|---|
| 179 | write(6,6) (x(k,j,1,3),k=1,lx1),(x(k,j,1,4),k=1,lx1)
|
|---|
| 180 | enddo
|
|---|
| 181 | write(6,*)
|
|---|
| 182 | write(6,*)
|
|---|
| 183 | endif
|
|---|
| 184 | c
|
|---|
| 185 | do j=ly1,1,-1
|
|---|
| 186 | write(6,6) (x(k,j,1,1),k=1,lx1),(x(k,j,1,2),k=1,lx1)
|
|---|
| 187 | enddo
|
|---|
| 188 | write(6,*)
|
|---|
| 189 | 3 format(4f6.2,5x,4f6.2)
|
|---|
| 190 | 6 format(6f8.4,5x,6f8.4)
|
|---|
| 191 | return
|
|---|
| 192 | end
|
|---|
| 193 | c-----------------------------------------------------------------------
|
|---|
| 194 | subroutine iarr_out(x,name)
|
|---|
| 195 | include 'SIZE'
|
|---|
| 196 | include 'INPUT'
|
|---|
| 197 | c
|
|---|
| 198 | integer x(lx1,ly1,lz1,lelt)
|
|---|
| 199 | character*13 name
|
|---|
| 200 | c
|
|---|
| 201 | if (nelv.gt.20) return
|
|---|
| 202 | write(6,*)
|
|---|
| 203 | write(6,1) name
|
|---|
| 204 | 1 format(a13)
|
|---|
| 205 | c
|
|---|
| 206 | c 3 D
|
|---|
| 207 | c
|
|---|
| 208 | if (if3d) then
|
|---|
| 209 | do iz=1,lz1
|
|---|
| 210 | write(6,*)
|
|---|
| 211 | do j=ly1,1,-1
|
|---|
| 212 | write(6,3) (x(k,j,iz,1),k=1,lx1),(x(k,j,iz,2),k=1,lx1)
|
|---|
| 213 | enddo
|
|---|
| 214 | enddo
|
|---|
| 215 | write(6,*)
|
|---|
| 216 | write(6,*)
|
|---|
| 217 | return
|
|---|
| 218 | endif
|
|---|
| 219 | c
|
|---|
| 220 | c 2 D
|
|---|
| 221 | c
|
|---|
| 222 | if (nelv.gt.2) then
|
|---|
| 223 | write(6,*)
|
|---|
| 224 | do j=ly1,1,-1
|
|---|
| 225 | write(6,6) (x(k,j,1,3),k=1,lx1),(x(k,j,1,4),k=1,lx1)
|
|---|
| 226 | enddo
|
|---|
| 227 | write(6,*)
|
|---|
| 228 | write(6,*)
|
|---|
| 229 | endif
|
|---|
| 230 | c
|
|---|
| 231 | do j=ly1,1,-1
|
|---|
| 232 | write(6,6) (x(k,j,1,1),k=1,lx1),(x(k,j,1,2),k=1,lx1)
|
|---|
| 233 | enddo
|
|---|
| 234 | write(6,*)
|
|---|
| 235 | 3 format(4i5,5x,4i5)
|
|---|
| 236 | 6 format(5i5,5x,5i5)
|
|---|
| 237 | return
|
|---|
| 238 | end
|
|---|
| 239 | c-----------------------------------------------------------------------
|
|---|
| 240 | subroutine iar2_out(x,name)
|
|---|
| 241 | include 'SIZE'
|
|---|
| 242 | c
|
|---|
| 243 | integer x(lx2,ly2,lz2,lelt)
|
|---|
| 244 | character*13 name
|
|---|
| 245 | c
|
|---|
| 246 | if (nelv.gt.20) return
|
|---|
| 247 | write(6,*)
|
|---|
| 248 | write(6,1) name
|
|---|
| 249 | 1 format(a13)
|
|---|
| 250 | if (nelv.gt.2) then
|
|---|
| 251 | write(6,*)
|
|---|
| 252 | do j=ly2,1,-1
|
|---|
| 253 | write(6,6) (x(k,j,1,3),k=1,lx2),(x(k,j,1,4),k=1,lx2)
|
|---|
| 254 | enddo
|
|---|
| 255 | write(6,*)
|
|---|
| 256 | write(6,*)
|
|---|
| 257 | endif
|
|---|
| 258 | c
|
|---|
| 259 | do j=ly2,1,-1
|
|---|
| 260 | write(6,6) (x(k,j,1,1),k=1,lx2),(x(k,j,1,2),k=1,lx2)
|
|---|
| 261 | enddo
|
|---|
| 262 | write(6,*)
|
|---|
| 263 | 6 format(3i5,5x,3i5)
|
|---|
| 264 | return
|
|---|
| 265 | end
|
|---|
| 266 | c-----------------------------------------------------------------------
|
|---|
| 267 | subroutine scsr_permute(bcsr,ib,jb,acsr,ia,ja,n
|
|---|
| 268 | $ ,icperm,inverse,nonzero,nzero)
|
|---|
| 269 | c
|
|---|
| 270 | c Permutes a Symmetric Compressed Sparse Row formatted matrix
|
|---|
| 271 | c to minimize bandwidth
|
|---|
| 272 | c
|
|---|
| 273 | c Reqired space: NONZERO is an array of size 2*nnz+2, where nnz
|
|---|
| 274 | c is the number of nonzeros. NONZERO and jb may
|
|---|
| 275 | c occupy the same space.
|
|---|
| 276 | c
|
|---|
| 277 | c NOTE!!: "inverse" is used as a scratch array for reals
|
|---|
| 278 | c in the swap call below... Therefore, inverse must
|
|---|
| 279 | c be aligned on even byte boundaries when running
|
|---|
| 280 | c in 64 bit precision.
|
|---|
| 281 | c
|
|---|
| 282 | c
|
|---|
| 283 | real bcsr(1)
|
|---|
| 284 | integer ib(1),jb(1)
|
|---|
| 285 | c
|
|---|
| 286 | real acsr(1)
|
|---|
| 287 | integer ia(1),ja(1)
|
|---|
| 288 | c
|
|---|
| 289 | integer icperm(n),nonzero(2,0:1)
|
|---|
| 290 | integer inverse(n),nzero(n)
|
|---|
| 291 | C HMT HACK ... if we're actualy using smbwr.c then we need to
|
|---|
| 292 | C uncomment the following line!!!
|
|---|
| 293 | C integer sm_bandwidth_reduction
|
|---|
| 294 | c
|
|---|
| 295 | integer icalld,max_n
|
|---|
| 296 | save icalld,max_n
|
|---|
| 297 | data icalld,max_n /0,0/
|
|---|
| 298 | C
|
|---|
| 299 | C HMT HACK ... if we're actualy using smbwr.c then we need to
|
|---|
| 300 | C delete the following line!!!
|
|---|
| 301 | write(6,*) 'HMT HACK in subroutine scsr_permute() ... pls fix!'
|
|---|
| 302 | call exitt
|
|---|
| 303 | C
|
|---|
| 304 | icalld=icalld+1
|
|---|
| 305 | c
|
|---|
| 306 | write(6,*) 'scsr_permute',n,acsr(1)
|
|---|
| 307 | mnz = 0
|
|---|
| 308 | mbw_in = 0
|
|---|
| 309 | do i=1,n
|
|---|
| 310 | n1 = ia(i) +1
|
|---|
| 311 | n2 = ia(i+1)-1
|
|---|
| 312 | do jj=n1,n2
|
|---|
| 313 | mnz = mnz+1
|
|---|
| 314 | j = ja(jj)
|
|---|
| 315 | nonzero(1,mnz) = i
|
|---|
| 316 | nonzero(2,mnz) = j
|
|---|
| 317 | C write (6,*) i,j
|
|---|
| 318 | mbw_in = max(mbw_in,abs(i-j))
|
|---|
| 319 | enddo
|
|---|
| 320 | enddo
|
|---|
| 321 | nonzero(1,0) = n
|
|---|
| 322 | nonzero(2,0) = mnz
|
|---|
| 323 | C
|
|---|
| 324 | itype = -1
|
|---|
| 325 | C itype = -2
|
|---|
| 326 | if (n.le.200) itype = -2
|
|---|
| 327 | write(6,*) 'this is mnz:',mnz
|
|---|
| 328 | C write(6,*) nonzero(1,0), nonzero(2,0)
|
|---|
| 329 | c
|
|---|
| 330 | C HMT HACK ... if we're actualy using smbwr.c then we need to
|
|---|
| 331 | C uncomment the following line!!!
|
|---|
| 332 | C ierr = sm_bandwidth_reduction(nonzero(1,0),icperm(1),itype)
|
|---|
| 333 | ierr = 1
|
|---|
| 334 | if (ierr.eq.0) then
|
|---|
| 335 | write(6,*) 'scsr_permute :: sm_bandwidth_reduction failed',ierr
|
|---|
| 336 | call exitt
|
|---|
| 337 | elseif (ierr.lt.0) then
|
|---|
| 338 | write(6,*) 'scsr_permute :: graph disconnected?',ierr
|
|---|
| 339 | call exitt
|
|---|
| 340 | endif
|
|---|
| 341 | c
|
|---|
| 342 | c Setup inverse map
|
|---|
| 343 | c
|
|---|
| 344 | do i=1,n
|
|---|
| 345 | inverse(icperm(i))=i
|
|---|
| 346 | enddo
|
|---|
| 347 | c
|
|---|
| 348 | C n20 = min(50,n)
|
|---|
| 349 | C write(6,20) (icperm (k),k=1,n20)
|
|---|
| 350 | C write(6,21) (inverse(k),k=1,n20)
|
|---|
| 351 | C 20 format('icp:',20i4)
|
|---|
| 352 | C 21 format('inv:',20i4)
|
|---|
| 353 | c
|
|---|
| 354 | c
|
|---|
| 355 | c Count *number* of nonzeros on each row in the new row pointer
|
|---|
| 356 | c
|
|---|
| 357 | call izero(nzero,2*mnz)
|
|---|
| 358 | do i=1,n
|
|---|
| 359 | n1 = ia(i)
|
|---|
| 360 | n2 = ia(i+1)-1
|
|---|
| 361 | ii = inverse(i)
|
|---|
| 362 | do jc=n1,n2
|
|---|
| 363 | j=ja(jc)
|
|---|
| 364 | jj = inverse(j)
|
|---|
| 365 | if (jj.ge.ii) then
|
|---|
| 366 | nzero(ii) = nzero(ii)+1
|
|---|
| 367 | else
|
|---|
| 368 | nzero(jj) = nzero(jj)+1
|
|---|
| 369 | endif
|
|---|
| 370 | enddo
|
|---|
| 371 | enddo
|
|---|
| 372 | c
|
|---|
| 373 | c Convert to pointers
|
|---|
| 374 | c
|
|---|
| 375 | ib(1) = 1
|
|---|
| 376 | do i=1,n
|
|---|
| 377 | ib(i+1) = ib(i)+nzero(i)
|
|---|
| 378 | enddo
|
|---|
| 379 | c
|
|---|
| 380 | c Fill permuted array
|
|---|
| 381 | c
|
|---|
| 382 | call icopy(nzero,ib,n)
|
|---|
| 383 | do i=1,n
|
|---|
| 384 | n1 = ia(i)
|
|---|
| 385 | n2 = ia(i+1)-1
|
|---|
| 386 | ii = inverse(i)
|
|---|
| 387 | do jc=n1,n2
|
|---|
| 388 | j=ja(jc)
|
|---|
| 389 | jj = inverse(j)
|
|---|
| 390 | if (jj.ge.ii) then
|
|---|
| 391 | jb (nzero(ii)) = jj
|
|---|
| 392 | bcsr(nzero(ii)) = acsr(jc)
|
|---|
| 393 | nzero(ii) = nzero(ii)+1
|
|---|
| 394 | else
|
|---|
| 395 | jb (nzero(jj)) = ii
|
|---|
| 396 | bcsr(nzero(jj)) = acsr(jc)
|
|---|
| 397 | nzero(jj) = nzero(jj)+1
|
|---|
| 398 | endif
|
|---|
| 399 | enddo
|
|---|
| 400 | enddo
|
|---|
| 401 | c
|
|---|
| 402 | c Sort each row of b in ascending column order --
|
|---|
| 403 | c just for subsequent access efficiency
|
|---|
| 404 | c
|
|---|
| 405 | m2 = mod(n2,2)+2
|
|---|
| 406 | do i=1,n
|
|---|
| 407 | n1 = ib(i)
|
|---|
| 408 | n2 = ib(i+1)-ib(i)
|
|---|
| 409 | call isort(jb (n1),nzero,n2)
|
|---|
| 410 | c write(6,*) 'call swap:',n,i,m2,n1,n2
|
|---|
| 411 | call swap (bcsr(n1),nzero,n2,inverse)
|
|---|
| 412 | enddo
|
|---|
| 413 | c
|
|---|
| 414 | if (icalld.le.10.or.mod(icalld,50).eq.0.or.n.gt.max_n) then
|
|---|
| 415 | mbw = mbw_csr(ib,jb,n)
|
|---|
| 416 | write(6,*) ierr,mbw,mbw_in,n,mnz,' New bandwidth'
|
|---|
| 417 | if (n.gt.max_n) max_n = n
|
|---|
| 418 | endif
|
|---|
| 419 | c
|
|---|
| 420 | return
|
|---|
| 421 | end
|
|---|
| 422 | c=======================================================================
|
|---|
| 423 | subroutine scsr_to_spb(abd,lda,acsr,ia,ja,n)
|
|---|
| 424 | c
|
|---|
| 425 | c This, and the companion routines map a *symmetrically*
|
|---|
| 426 | c stored (upper 1/2 only) Compressed Sparse Row formatted
|
|---|
| 427 | c matrix to the linpack banded format.
|
|---|
| 428 | c
|
|---|
| 429 | c
|
|---|
| 430 | real abd(1),acsr(1)
|
|---|
| 431 | integer m,n,ia(1),ja(1)
|
|---|
| 432 | c
|
|---|
| 433 | lda = mbw_csr(ia,ja,n) + 1
|
|---|
| 434 | call scsr_to_spbm(abd,lda,acsr,ia,ja,n)
|
|---|
| 435 | c
|
|---|
| 436 | return
|
|---|
| 437 | end
|
|---|
| 438 | c=======================================================================
|
|---|
| 439 | subroutine scsr_to_spbm(abd,lda,acsr,ia,ja,n)
|
|---|
| 440 | c
|
|---|
| 441 | real abd(lda,1),acsr(1)
|
|---|
| 442 | integer n,ia(1),ja(1)
|
|---|
| 443 | c
|
|---|
| 444 | call rzero(abd,lda*n)
|
|---|
| 445 | do i=1,n
|
|---|
| 446 | n1 = ia(i)
|
|---|
| 447 | n2 = ia(i+1)-1
|
|---|
| 448 | do jc=n1,n2
|
|---|
| 449 | j=ja(jc)
|
|---|
| 450 | klin = lda + (i-j)
|
|---|
| 451 | jlin = j
|
|---|
| 452 | abd(klin,jlin) = acsr(jc)
|
|---|
| 453 | enddo
|
|---|
| 454 | enddo
|
|---|
| 455 | return
|
|---|
| 456 | end
|
|---|
| 457 | function mbw_csr(ia,ja,n)
|
|---|
| 458 | c
|
|---|
| 459 | c Determine the maximum bandwidth for a csr formatted matrix
|
|---|
| 460 | c
|
|---|
| 461 | integer ia(1),ja(1)
|
|---|
| 462 | c
|
|---|
| 463 | m = 0
|
|---|
| 464 | do i=1,n
|
|---|
| 465 | n1 = ia(i)
|
|---|
| 466 | n2 = ia(i+1)-1
|
|---|
| 467 | do jc=n1,n2
|
|---|
| 468 | j=ja(jc)
|
|---|
| 469 | m = max(m,abs(i-j))
|
|---|
| 470 | enddo
|
|---|
| 471 | enddo
|
|---|
| 472 | mbw_csr = m
|
|---|
| 473 | return
|
|---|
| 474 | end
|
|---|
| 475 | c=======================================================================
|
|---|
| 476 | subroutine out_spbmat(abd,n,lda,name)
|
|---|
| 477 | real abd(lda,1)
|
|---|
| 478 | c
|
|---|
| 479 | character*9 name
|
|---|
| 480 | character*6 s(22)
|
|---|
| 481 | c
|
|---|
| 482 | m = lda-1
|
|---|
| 483 | write(6,1) name,n,m
|
|---|
| 484 | 1 format(/,'SPB Mat:',a9,3x,'n =',i3,' m =',i3,/)
|
|---|
| 485 | c
|
|---|
| 486 | n22 = min(n,22)
|
|---|
| 487 | do i=1,n
|
|---|
| 488 | call blank(s,132)
|
|---|
| 489 | ii = lda*i
|
|---|
| 490 | do k=0,m
|
|---|
| 491 | j = i+k
|
|---|
| 492 | a = abd(ii+k*m,1)
|
|---|
| 493 | if (a.ne.0..and.j.le.n22) write(s(j),6) a
|
|---|
| 494 | enddo
|
|---|
| 495 | write(6,22) (s(k),k=1,n22)
|
|---|
| 496 | enddo
|
|---|
| 497 | 6 format(f6.3)
|
|---|
| 498 | 22 format(22a6)
|
|---|
| 499 | c
|
|---|
| 500 | return
|
|---|
| 501 | end
|
|---|
| 502 | c=======================================================================
|
|---|
| 503 | subroutine swap(b,ind,n,temp)
|
|---|
| 504 | real B(1),TEMP(1)
|
|---|
| 505 | integer IND(1)
|
|---|
| 506 | C***
|
|---|
| 507 | C*** SORT ASSOCIATED ELEMENTS BY PUTTING ITEM(JJ)
|
|---|
| 508 | C*** INTO ITEM(I), WHERE JJ=IND(I).
|
|---|
| 509 | C***
|
|---|
| 510 | DO 20 I=1,N
|
|---|
| 511 | JJ=IND(I)
|
|---|
| 512 | TEMP(I)=B(JJ)
|
|---|
| 513 | 20 CONTINUE
|
|---|
| 514 | DO 30 I=1,N
|
|---|
| 515 | 30 B(I)=TEMP(I)
|
|---|
| 516 | RETURN
|
|---|
| 517 | END
|
|---|
| 518 | c=======================================================================
|
|---|
| 519 | subroutine ipermute(a,icperm,n,b)
|
|---|
| 520 | integer a(1),b(1)
|
|---|
| 521 | integer icperm(1)
|
|---|
| 522 | c
|
|---|
| 523 | call icopy(b,a,n)
|
|---|
| 524 | do i=1,n
|
|---|
| 525 | a(i) = b(icperm(i))
|
|---|
| 526 | enddo
|
|---|
| 527 | return
|
|---|
| 528 | end
|
|---|
| 529 | c=======================================================================
|
|---|
| 530 | subroutine out_csrmat(acsr,ia,ja,n,name9)
|
|---|
| 531 | real acsr(1)
|
|---|
| 532 | integer ia(1),ja(1)
|
|---|
| 533 | c
|
|---|
| 534 | character*9 name9
|
|---|
| 535 | character*6 s(22)
|
|---|
| 536 | c
|
|---|
| 537 | nnz = ia(n+1)-ia(1)
|
|---|
| 538 | c
|
|---|
| 539 | write(6,1) name9,n,nnz
|
|---|
| 540 | 1 format(/,'CSR Mat:',a9,3x,'n =',i3,3x,'nnz =',i5,/)
|
|---|
| 541 | c
|
|---|
| 542 | n22 = min(n,22)
|
|---|
| 543 | n29 = min(n,29)
|
|---|
| 544 | do i=1,n29
|
|---|
| 545 | call blank(s,132)
|
|---|
| 546 | n1 = ia(i)
|
|---|
| 547 | n2 = ia(i+1)-1
|
|---|
| 548 | do jj=n1,n2
|
|---|
| 549 | j = ja (jj)
|
|---|
| 550 | a = acsr(jj)
|
|---|
| 551 | if (a.ne.0..and.j.le.n22) write(s(j),6) a
|
|---|
| 552 | enddo
|
|---|
| 553 | write(6,22) (s(k),k=1,n22)
|
|---|
| 554 | enddo
|
|---|
| 555 | 6 format(f6.2)
|
|---|
| 556 | 22 format(22a6)
|
|---|
| 557 | c
|
|---|
| 558 | return
|
|---|
| 559 | end
|
|---|
| 560 | c=======================================================================
|
|---|
