EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
genmat.f
Go to the documentation of this file.
1 c----------------------------------------------------------------------c
2 c S P A R S K I T c
3 c----------------------------------------------------------------------c
4 c MATRIX GENERATION ROUTINES -- FINITE DIFFERENCE MATRICES c
5 c----------------------------------------------------------------------c
6 c contents: c
7 c---------- c
8 c gen57pt : generates 5-point and 7-point matrices. c
9 c gen57bl : generates block 5-point and 7-point matrices. c
10 c c
11 c supporting routines: c
12 c--------- c
13 c gensten : generate the stencil (point version) c
14 c bsten : generate the stencil (block version) c
15 c fdaddbc : finite difference add boundary conditions c
16 c fdreduce : reduce the system to eliminate node with known values c
17 c clrow : clear a row of a CSR matrix c
18 c lctcsr : locate the position of A(i,j) in CSR format c
19 c----------------------------------------------------------------------c
20  subroutine gen57pt(nx,ny,nz,al,mode,n,a,ja,ia,iau,rhs)
21  integer ja(*),ia(*),iau(*), nx, ny, nz, mode, n
22  real*8 a(*), rhs(*), al(6)
23 c-----------------------------------------------------------------------
24 c On entry:
25 c
26 c nx = number of grid points in x direction
27 c ny = number of grid points in y direction
28 c nz = number of grid points in z direction
29 c al = array of size 6, carries the coefficient alpha of the
30 c boundary conditions
31 c mode = what to generate:
32 c < 0 : generate the graph only,
33 c = 0 : generate the matrix,
34 c > 0 : generate the matrix and the right-hand side.
35 c
36 c On exit:
37 c
38 c n = number of nodes with unknown values, ie number of rows
39 c in the matrix
40 c
41 c a,ja,ia = resulting matrix in row-sparse format
42 c
43 c iau = integer*n, containing the poisition of the diagonal element
44 c in the a, ja, ia structure
45 c
46 c rhs = the right-hand side
47 c
48 c External functions needed (must be supplied by caller)
49 c afun, bfun, cfun, dfun, efun, ffun, gfun, hfun
50 c betfun, gamfun
51 c They have the following prototype:
52 c real*8 function xfun(x, y, z)
53 c real*8 x, y, z
54 c-----------------------------------------------------------------------
55 c This subroutine computes the sparse matrix in compressed sparse row
56 c format for the elliptic equation:
57 c d du d du d du du du du
58 c L u = --(A --) + --(B --) + --(C --) + D -- + E -- + F -- + G u = H u
59 c dx dx dy dy dz dz dx dy dz
60 c
61 c with general Mixed Boundary conditions, on a rectangular 1-D,
62 c 2-D or 3-D grid using 2nd order centered difference schemes.
63 c
64 c The functions a, b, ..., g, h are known through the
65 c as afun, bfun, ..., gfun, hfun in this subroutine.
66 c NOTE: To obtain the correct matrix, any function that is not
67 c needed should be set to zero. For example for two-dimensional
68 c problems, nz should be set to 1 and the functions cfun and ffun
69 c should be zero functions.
70 c
71 c The Boundary condition is specified in the following form:
72 c du
73 c alpha -- + beta u = gamma
74 c dn
75 c Where alpha is constant at each side of the boundary surfaces. Alpha
76 c is represented by parameter al. It is expected to an array that
77 c contains enough elements to specify the boundaries for the problem,
78 c 1-D case needs two elements, 2-D needs 4 and 3-D needs 6. The order
79 c of the boundaries in the array is left(west), right(east),
80 c bottom(south), top(north), front, rear. Beta and gamma are functions
81 c of type real with three arguments x, y, z. These two functions are
82 c known subroutine 'addbc' as betfun and gamfun. They should following
83 c the same notion as afun ... hfun. For more restriction on afun ...
84 c hfun, please read the documentation follows the subroutine 'getsten',
85 c and, for more on betfun and gamfun, please refer to the documentation
86 c under subroutine 'fdaddbc'.
87 c
88 c The nodes are ordered using natural ordering, first x direction, then
89 c y, then z. The mesh size h is uniform and determined by grid points
90 c in the x-direction.
91 c
92 c The domain specified for the problem is [0 .ge. x .ge. 1],
93 c [0 .ge. y .ge. (ny-1)*h] and [0 .ge. z .ge. (nz-1)*h], where h is
94 c 1 / (nx-1). Thus if non-Dirichlet boundary condition is specified,
95 c the mesh will have nx points along the x direction, ny along y and
96 c nz along z. For 1-D case, both y and z value are assumed to zero
97 c when calling relavent functions that have three parameters.
98 c Similarly, for 2-D case, z is assumed to be zero.
99 c
100 c About the expectation of nx, ny and nz:
101 c nx is required to be .gt. 1 always;
102 c if the second dimension is present in the problem, then ny should be
103 c .gt. 1, else 1;
104 c if the third dimension is present in the problem, nz .gt. 1, else 1.
105 c when ny is 1, nz must be 1.
106 c-----------------------------------------------------------------------
107 c
108 c stencil [1:7] has the following meaning:
109 c
110 c center point = stencil(1)
111 c west point = stencil(2)
112 c east point = stencil(3)
113 c south point = stencil(4)
114 c north point = stencil(5)
115 c front point = stencil(6)
116 c back point = stencil(7)
117 c
118 c al[1:6] carry the coefficient alpha in the similar order
119 c
120 c west side = al(1)
121 c east side = al(2)
122 c south side = al(3)
123 c north side = al(4)
124 c front side = al(5)
125 c back side = al(6)
126 c
127 c al(4)
128 c st(5)
129 c |
130 c |
131 c | al(6)
132 c | .st(7)
133 c | .
134 c al(1) | . al(2)
135 c st(2) ----------- st(1) ---------- st(3)
136 c . |
137 c . |
138 c . |
139 c st(6) |
140 c al(5) |
141 c |
142 c st(4)
143 c al(3)
144 c
145 c-------------------------------------------------------------------
146 c some constants
147 c
148  real*8 one
149  parameter(one=1.0d0)
150 c
151 c local variables
152 c
153  integer ix, iy, iz, kx, ky, kz, node, iedge
154  real*8 r, h, stencil(7)
155  logical value, genrhs
156 c
157 c nx has to be larger than 1
158 c
159  if (nx.le.1) return
160  h = one / dble(nx-1)
161 c
162 c the mode
163 c
164  value = (mode.ge.0)
165  genrhs = (mode.gt.0)
166 c
167 c first generate the whole matrix as if the boundary condition does
168 c not exist
169 c
170  kx = 1
171  ky = nx
172  kz = nx*ny
173  iedge = 1
174  node = 1
175  do 100 iz = 1,nz
176  do 90 iy = 1,ny
177  do 80 ix = 1,nx
178  ia(node) = iedge
179 c
180 c compute the stencil at the current node
181 c
182  if (value) call
183  & getsten(nx,ny,nz,mode,ix-1,iy-1,iz-1,stencil,h,r)
184 c west
185  if (ix.gt.1) then
186  ja(iedge)=node-kx
187  if (value) a(iedge) = stencil(2)
188  iedge=iedge + 1
189  end if
190 c south
191  if (iy.gt.1) then
192  ja(iedge)=node-ky
193  if (value) a(iedge) = stencil(4)
194  iedge=iedge + 1
195  end if
196 c front plane
197  if (iz.gt.1) then
198  ja(iedge)=node-kz
199  if (value) a(iedge) = stencil(6)
200  iedge=iedge + 1
201  endif
202 c center node
203  ja(iedge) = node
204  iau(node) = iedge
205  if (value) a(iedge) = stencil(1)
206  iedge = iedge + 1
207 c east
208  if (ix.lt.nx) then
209  ja(iedge)=node+kx
210  if (value) a(iedge) = stencil(3)
211  iedge=iedge + 1
212  end if
213 c north
214  if (iy.lt.ny) then
215  ja(iedge)=node+ky
216  if (value) a(iedge) = stencil(5)
217  iedge=iedge + 1
218  end if
219 c back plane
220  if (iz.lt.nz) then
221  ja(iedge)=node+kz
222  if (value) a(iedge) = stencil(7)
223  iedge=iedge + 1
224  end if
225 c the right-hand side
226  if (genrhs) rhs(node) = r
227  node=node+1
228  80 continue
229  90 continue
230  100 continue
231  ia(node)=iedge
232 c
233 c Add in the boundary conditions
234 c
235  call fdaddbc(nx,ny,nz,a,ja,ia,iau,rhs,al,h)
236 c
237 c eliminate the boudary nodes from the matrix
238 c
239  call fdreduce(nx,ny,nz,al,n,a,ja,ia,iau,rhs,stencil)
240 c
241 c done
242 c
243  return
244 c-----end-of-gen57pt----------------------------------------------------
245 c-----------------------------------------------------------------------
246  end
247 c-----------------------------------------------------------------------
248  subroutine getsten (nx,ny,nz,mode,kx,ky,kz,stencil,h,rhs)
249  integer nx,ny,nz,mode,kx,ky,kz
250  real*8 stencil(*),h,rhs,afun,bfun,cfun,dfun,efun,ffun,gfun,hfun
251  external afun,bfun,cfun,dfun,efun,ffun,gfun,hfun
252 c-----------------------------------------------------------------------
253 c This subroutine calculates the correct stencil values for
254 c centered difference discretization of the elliptic operator
255 c and the right-hand side
256 c
257 c L u = delx( A delx u ) + dely ( B dely u) + delz ( C delz u ) +
258 c delx ( D u ) + dely (E u) + delz( F u ) + G u = H
259 c
260 c For 2-D problems the discretization formula that is used is:
261 c
262 c h**2 * Lu == A(i+1/2,j)*{u(i+1,j) - u(i,j)} +
263 c A(i-1/2,j)*{u(i-1,j) - u(i,j)} +
264 c B(i,j+1/2)*{u(i,j+1) - u(i,j)} +
265 c B(i,j-1/2)*{u(i,j-1) - u(i,j)} +
266 c (h/2)*D(i,j)*{u(i+1,j) - u(i-1,j)} +
267 c (h/2)*E(i,j)*{u(i,j+1) - u(i,j-1)} +
268 c (h/2)*E(i,j)*{u(i,j+1) - u(i,j-1)} +
269 c (h**2)*G(i,j)*u(i,j)
270 c-----------------------------------------------------------------------
271 c some constants
272 c
273  real*8 zero, half
274  parameter(zero=0.0d0,half=0.5d0)
275 c
276 c local variables
277 c
278  integer k
279  real*8 hhalf,cntr, x, y, z, coeff
280 c
281 c if mode < 0, we shouldn't have come here
282 c
283  if (mode .lt. 0) return
284 c
285  do 200 k=1,7
286  stencil(k) = zero
287  200 continue
288 c
289  hhalf = h*half
290  x = h*dble(kx)
291  y = h*dble(ky)
292  z = h*dble(kz)
293  cntr = zero
294 c differentiation wrt x:
295  coeff = afun(x+hhalf,y,z)
296  stencil(3) = stencil(3) + coeff
297  cntr = cntr + coeff
298 c
299  coeff = afun(x-hhalf,y,z)
300  stencil(2) = stencil(2) + coeff
301  cntr = cntr + coeff
302 c
303  coeff = dfun(x,y,z)*hhalf
304  stencil(3) = stencil(3) + coeff
305  stencil(2) = stencil(2) - coeff
306  if (ny .le. 1) goto 99
307 c
308 c differentiation wrt y:
309 c
310  coeff = bfun(x,y+hhalf,z)
311  stencil(5) = stencil(5) + coeff
312  cntr = cntr + coeff
313 c
314  coeff = bfun(x,y-hhalf,z)
315  stencil(4) = stencil(4) + coeff
316  cntr = cntr + coeff
317 c
318  coeff = efun(x,y,z)*hhalf
319  stencil(5) = stencil(5) + coeff
320  stencil(4) = stencil(4) - coeff
321  if (nz .le. 1) goto 99
322 c
323 c differentiation wrt z:
324 c
325  coeff = cfun(x,y,z+hhalf)
326  stencil(7) = stencil(7) + coeff
327  cntr = cntr + coeff
328 c
329  coeff = cfun(x,y,z-hhalf)
330  stencil(6) = stencil(6) + coeff
331  cntr = cntr + coeff
332 c
333  coeff = ffun(x,y,z)*hhalf
334  stencil(7) = stencil(7) + coeff
335  stencil(6) = stencil(6) - coeff
336 c
337 c contribution from function G:
338 c
339  99 coeff = gfun(x,y,z)
340  stencil(1) = h*h*coeff - cntr
341 c
342 c the right-hand side
343 c
344  if (mode .gt. 0) rhs = h*h*hfun(x,y,z)
345 c
346  return
347 c------end-of-getsten---------------------------------------------------
348 c-----------------------------------------------------------------------
349  end
350 c-----------------------------------------------------------------------
351  subroutine gen57bl (nx,ny,nz,nfree,na,n,a,ja,ia,iau,stencil)
352 c implicit real*8 (a-h,o-z)
353  integer ja(*),ia(*),iau(*),nx,ny,nz,nfree,na,n
354  real*8 a(na,1), stencil(7,1)
355 c--------------------------------------------------------------------
356 c This subroutine computes the sparse matrix in compressed
357 c format for the elliptic operator
358 c
359 c L u = delx( a . delx u ) + dely ( b . dely u) + delz ( c . delz u ) +
360 c delx ( d . u ) + dely (e . u) + delz( f . u ) + g . u
361 c
362 c Here u is a vector of nfree componebts and each of the functions
363 c a, b, c, d, e, f, g is an (nfree x nfree) matrix depending of
364 c the coordinate (x,y,z).
365 c with Dirichlet Boundary conditions, on a rectangular 1-D,
366 c 2-D or 3-D grid using centered difference schemes.
367 c
368 c The functions a, b, ..., g are known through the
369 c subroutines afunbl, bfunbl, ..., gfunbl. (user supplied) .
370 c
371 c uses natural ordering, first x direction, then y, then z
372 c mesh size h is uniform and determined by grid points
373 c in the x-direction.
374 c
375 c The output matrix is in Block -- Sparse Row format.
376 c
377 c--------------------------------------------------------------------
378 c parameters:
379 c-------------
380 c Input:
381 c ------
382 c nx = number of points in x direction
383 c ny = number of points in y direction
384 c nz = number of points in z direction
385 c nfree = number of degrees of freedom per point
386 c na = first dimension of array a as declared in calling
387 c program. Must be .ge. nfree**2
388 c
389 c Output:
390 c ------
391 c n = dimension of matrix (output)
392 c
393 c a, ja, ia = resulting matrix in Block Sparse Row format
394 c a(1:nfree**2, j ) contains a nonzero block and ja(j)
395 c contains the (block) column number of this block.
396 c the block dimension of the matrix is n (output) and
397 c therefore the total number of (scalar) rows is n x nfree.
398 c
399 c iau = integer*n containing the position of the diagonal element
400 c in the a, ja, ia structure
401 c
402 c Work space:
403 c------------
404 c stencil = work array of size (7,nfree**2) [stores local stencils]
405 c
406 c--------------------------------------------------------------------
407 c
408 c stencil (1:7,*) has the following meaning:
409 c
410 c center point = stencil(1)
411 c west point = stencil(2)
412 c east point = stencil(3)
413 c south point = stencil(4)
414 c north point = stencil(5)
415 c front point = stencil(6)
416 c back point = stencil(7)
417 c
418 c
419 c st(5)
420 c |
421 c |
422 c |
423 c | .st(7)
424 c | .
425 c | .
426 c st(2) ----------- st(1) ---------- st(3)
427 c . |
428 c . |
429 c . |
430 c st(6) |
431 c |
432 c |
433 c st(4)
434 c
435 c-------------------------------------------------------------------
436 c some constants
437 c
438  real*8 one
439  parameter(one=1.0d0)
440 c
441 c local variables
442 c
443  integer iedge,ix,iy,iz,k,kx,ky,kz,nfree2,node
444  real*8 h
445 c
446  h = one/dble(nx+1)
447  kx = 1
448  ky = nx
449  kz = nx*ny
450  nfree2 = nfree*nfree
451  iedge = 1
452  node = 1
453  do 100 iz = 1,nz
454  do 90 iy = 1,ny
455  do 80 ix = 1,nx
456  ia(node) = iedge
457  call bsten(nx,ny,nz,ix,iy,iz,nfree,stencil,h)
458 c west
459  if (ix.gt.1) then
460  ja(iedge)=node-kx
461  do 4 k=1,nfree2
462  a(k,iedge) = stencil(2,k)
463  4 continue
464  iedge=iedge + 1
465  end if
466 c south
467  if (iy.gt.1) then
468  ja(iedge)=node-ky
469  do 5 k=1,nfree2
470  a(k,iedge) = stencil(4,k)
471  5 continue
472  iedge=iedge + 1
473  end if
474 c front plane
475  if (iz.gt.1) then
476  ja(iedge)=node-kz
477  do 6 k=1,nfree2
478  a(k,iedge) = stencil(6,k)
479  6 continue
480  iedge=iedge + 1
481  endif
482 c center node
483  ja(iedge) = node
484  iau(node) = iedge
485  do 7 k=1,nfree2
486  a(k,iedge) = stencil(1,k)
487  7 continue
488  iedge = iedge + 1
489 c -- upper part
490 c east
491  if (ix.lt.nx) then
492  ja(iedge)=node+kx
493  do 8 k=1,nfree2
494  a(k,iedge) = stencil(3,k)
495  8 continue
496  iedge=iedge + 1
497  end if
498 c north
499  if (iy.lt.ny) then
500  ja(iedge)=node+ky
501  do 9 k=1,nfree2
502  a(k,iedge) = stencil(5,k)
503  9 continue
504  iedge=iedge + 1
505  end if
506 c back plane
507  if (iz.lt.nz) then
508  ja(iedge)=node+kz
509  do 10 k=1,nfree2
510  a(k,iedge) = stencil(7,k)
511  10 continue
512  iedge=iedge + 1
513  end if
514 c------next node -------------------------
515  node=node+1
516  80 continue
517  90 continue
518  100 continue
519 c
520 c -- new version of BSR -- renumbering removed.
521 c change numbering of nodes so that each ja(k) will contain the
522 c actual column number in the original matrix of entry (1,1) of each
523 c block (k).
524 c do 101 k=1,iedge-1
525 c ja(k) = (ja(k)-1)*nfree+1
526 c 101 continue
527 c
528 c n = (node-1)*nfree
529  n = node-1
530  ia(node)=iedge
531  return
532 c--------------end-of-gen57bl-------------------------------------------
533 c-----------------------------------------------------------------------
534  end
535 c-----------------------------------------------------------------------
536  subroutine bsten (nx,ny,nz,kx,ky,kz,nfree,stencil,h)
537 c-----------------------------------------------------------------------
538 c This subroutine calcultes the correct block-stencil values for
539 c centered difference discretization of the elliptic operator
540 c (block version of stencil)
541 c
542 c L u = delx( a delx u ) + dely ( b dely u) + delz ( c delz u ) +
543 c d delx ( u ) + e dely (u) + f delz( u ) + g u
544 c
545 c For 2-D problems the discretization formula that is used is:
546 c
547 c h**2 * Lu == a(i+1/2,j)*{u(i+1,j) - u(i,j)} +
548 c a(i-1/2,j)*{u(i-1,j) - u(i,j)} +
549 c b(i,j+1/2)*{u(i,j+1) - u(i,j)} +
550 c b(i,j-1/2)*{u(i,j-1) - u(i,j)} +
551 c (h/2)*d(i,j)*{u(i+1,j) - u(i-1,j)} +
552 c (h/2)*e(i,j)*{u(i,j+1) - u(i,j-1)} +
553 c (h/2)*e(i,j)*{u(i,j+1) - u(i,j-1)} +
554 c (h**2)*g(i,j)*u(i,j)
555 c-----------------------------------------------------------------------
556 c some constants
557 c
558  real*8 zero,half
559  parameter(zero=0.0d0,half=0.5d0)
560 c
561 c local variables
562 c
563  integer i,k,kx,ky,kz,nfree,nfree2,nx,ny,nz
564  real*8 stencil(7,*), cntr(225), coeff(225),h,h2,hhalf,x,y,z
565 c------------
566  if (nfree .gt. 15) then
567  print *, ' ERROR ** nfree too large '
568  stop
569  endif
570 c
571  nfree2 = nfree*nfree
572  do 200 k=1, nfree2
573  cntr(k) = zero
574  do 199 i=1,7
575  stencil(i,k) = zero
576  199 continue
577  200 continue
578 c------------
579  hhalf = h*half
580  h2 = h*h
581  x = h*dble(kx)
582  y = h*dble(ky)
583  z = h*dble(kz)
584 c differentiation wrt x:
585  call afunbl(nfree,x+hhalf,y,z,coeff)
586  do 1 k=1, nfree2
587  stencil(3,k) = stencil(3,k) + coeff(k)
588  cntr(k) = cntr(k) + coeff(k)
589  1 continue
590 c
591  call afunbl(nfree,x-hhalf,y,z,coeff)
592  do 2 k=1, nfree2
593  stencil(2,k) = stencil(2,k) + coeff(k)
594  cntr(k) = cntr(k) + coeff(k)
595  2 continue
596 c
597  call dfunbl(nfree,x,y,z,coeff)
598  do 3 k=1, nfree2
599  stencil(3,k) = stencil(3,k) + coeff(k)*hhalf
600  stencil(2,k) = stencil(2,k) - coeff(k)*hhalf
601  3 continue
602  if (ny .le. 1) goto 99
603 c
604 c differentiation wrt y:
605 c
606  call bfunbl(nfree,x,y+hhalf,z,coeff)
607  do 4 k=1,nfree2
608  stencil(5,k) = stencil(5,k) + coeff(k)
609  cntr(k) = cntr(k) + coeff(k)
610  4 continue
611 c
612  call bfunbl(nfree,x,y-hhalf,z,coeff)
613  do 5 k=1, nfree2
614  stencil(4,k) = stencil(4,k) + coeff(k)
615  cntr(k) = cntr(k) + coeff(k)
616  5 continue
617 c
618  call efunbl(nfree,x,y,z,coeff)
619  do 6 k=1, nfree2
620  stencil(5,k) = stencil(5,k) + coeff(k)*hhalf
621  stencil(4,k) = stencil(4,k) - coeff(k)*hhalf
622  6 continue
623  if (nz .le. 1) goto 99
624 c
625 c differentiation wrt z:
626 c
627  call cfunbl(nfree,x,y,z+hhalf,coeff)
628  do 7 k=1, nfree2
629  stencil(7,k) = stencil(7,k) + coeff(k)
630  cntr(k) = cntr(k) + coeff(k)
631  7 continue
632 c
633  call cfunbl(nfree,x,y,z-hhalf,coeff)
634  do 8 k=1, nfree2
635  stencil(6,k) = stencil(6,k) + coeff(k)
636  cntr(k) = cntr(k) + coeff(k)
637  8 continue
638 c
639  call ffunbl(nfree,x,y,z,coeff)
640  do 9 k=1, nfree2
641  stencil(7,k) = stencil(7,k) + coeff(k)*hhalf
642  stencil(6,k) = stencil(6,k) - coeff(k)*hhalf
643  9 continue
644 c
645 c discretization of product by g:
646 c
647  99 call gfunbl(nfree,x,y,z,coeff)
648  do 10 k=1, nfree2
649  stencil(1,k) = h2*coeff(k) - cntr(k)
650  10 continue
651 c
652  return
653 c------------end of bsten-----------------------------------------------
654 c-----------------------------------------------------------------------
655  end
656  subroutine fdreduce(nx,ny,nz,alpha,n,a,ja,ia,iau,rhs,stencil)
657  implicit none
658  integer nx,ny, nz, n, ia(*), ja(*), iau(*)
659  real*8 alpha(*), a(*), rhs(*), stencil(*)
660 c-----------------------------------------------------------------------
661 c This subroutine tries to reduce the size of the matrix by looking
662 c for Dirichlet boundary conditions at each surface and solve the boundary
663 c value and modify the right-hand side of related nodes, then clapse all
664 c the boundary nodes.
665 c-----------------------------------------------------------------------
666 c parameters
667 c
668  real*8 zero
669  parameter(zero=0.0d0)
670 c
671 c local variables
672 c
673  integer i,j,k,kx,ky,kz,lx,ux,ly,uy,lz,uz,node,nbnode,lk,ld,iedge
674  real*8 val
675  integer lctcsr
676  external lctcsr
677 c
678 c The first half of this subroutine will try to change the right-hand
679 c side of all the nodes that has a neighbor with Dirichlet boundary
680 c condition, since in this case the value of the boundary point is
681 c known.
682 c Then in the second half, we will try to eliminate the boundary
683 c points with known values (with Dirichlet boundary condition).
684 c
685  kx = 1
686  ky = nx
687  kz = nx*ny
688  lx = 1
689  ux = nx
690  ly = 1
691  uy = ny
692  lz = 1
693  uz = nz
694 c
695 c Here goes the first part. ----------------------------------------
696 c
697 c the left (west) side
698 c
699  if (alpha(1) .eq. zero) then
700  lx = 2
701  do 10 k = 1, nz
702  do 11 j = 1, ny
703  node = (k-1)*kz + (j-1)*ky + 1
704  nbnode = node + kx
705  lk = lctcsr(nbnode, node, ja, ia)
706  ld = iau(node)
707  val = rhs(node)/a(ld)
708 c modify the rhs
709  rhs(nbnode) = rhs(nbnode) - a(lk)*val
710  11 continue
711  10 continue
712  endif
713 c
714 c right (east) side
715 c
716  if (alpha(2) .eq. zero) then
717  ux = nx - 1
718  do 20 k = 1, nz
719  do 21 j = 1, ny
720  node = (k-1)*kz + (j-1)*ky + nx
721  nbnode = node - kx
722  lk = lctcsr(nbnode, node, ja, ia)
723  ld = iau(node)
724  val = rhs(node)/a(ld)
725 c modify the rhs
726  rhs(nbnode) = rhs(nbnode) - a(lk)*val
727  21 continue
728  20 continue
729  endif
730 c
731 c if it's only 1-D, skip the following part
732 c
733  if (ny .le. 1) goto 100
734 c
735 c the bottom (south) side
736 c
737  if (alpha(3) .eq. zero) then
738  ly = 2
739  do 30 k = 1, nz
740  do 31 i = lx, ux
741  node = (k-1)*kz + i
742  nbnode = node + ky
743  lk = lctcsr(nbnode, node, ja, ia)
744  ld = iau(node)
745  val = rhs(node)/a(ld)
746 c modify the rhs
747  rhs(nbnode) = rhs(nbnode) - a(lk)*val
748  31 continue
749  30 continue
750  endif
751 c
752 c top (north) side
753 c
754  if (alpha(4) .eq. zero) then
755  uy = ny - 1
756  do 40 k = 1, nz
757  do 41 i = lx, ux
758  node = (k-1)*kz + i + (ny-1)*ky
759  nbnode = node - ky
760  lk = lctcsr(nbnode, node, ja, ia)
761  ld = iau(node)
762  val = rhs(node)/a(ld)
763 c modify the rhs
764  rhs(nbnode) = rhs(nbnode) - a(lk)*val
765  41 continue
766  40 continue
767  endif
768 c
769 c if only 2-D skip the following section on z
770 c
771  if (nz .le. 1) goto 100
772 c
773 c the front surface
774 c
775  if (alpha(5) .eq. zero) then
776  lz = 2
777  do 50 j = ly, uy
778  do 51 i = lx, ux
779  node = (j-1)*ky + i
780  nbnode = node + kz
781  lk = lctcsr(nbnode, node, ja, ia)
782  ld = iau(node)
783  val = rhs(node)/a(ld)
784 c modify the rhs
785  rhs(nbnode) = rhs(nbnode) - a(lk)*val
786  51 continue
787  50 continue
788  endif
789 c
790 c rear surface
791 c
792  if (alpha(6) .eq. zero) then
793  uz = nz - 1
794  do 60 j = ly, uy
795  do 61 i = lx, ux
796  node = (nz-1)*kz + (j-1)*ky + i
797  nbnode = node - kz
798  lk = lctcsr(nbnode, node, ja, ia)
799  ld = iau(node)
800  val = rhs(node)/a(ld)
801 c modify the rhs
802  rhs(nbnode) = rhs(nbnode) - a(lk)*val
803  61 continue
804  60 continue
805  endif
806 c
807 c now the second part ----------------------------------------------
808 c
809 c go through all the actual nodes with unknown values, collect all
810 c of them to form a new matrix in compressed sparse row format.
811 c
812  100 kx = 1
813  ky = ux - lx + 1
814  kz = (uy - ly + 1) * ky
815  node = 1
816  iedge = 1
817  do 80 k = lz, uz
818  do 81 j = ly, uy
819  do 82 i = lx, ux
820 c
821 c the corresponding old node number
822  nbnode = ((k-1)*ny + j-1)*nx + i
823 c
824 c copy the row into local stencil, copy is done is the exact
825 c same order as the stencil is written into array a
826  lk = ia(nbnode)
827  if (i.gt.1) then
828  stencil(2) = a(lk)
829  lk = lk + 1
830  end if
831  if (j.gt.1) then
832  stencil(4) = a(lk)
833  lk = lk + 1
834  end if
835  if (k.gt.1) then
836  stencil(6) = a(lk)
837  lk = lk + 1
838  end if
839  stencil(1) = a(lk)
840  lk = lk + 1
841  if (i.lt.nx) then
842  stencil(3) = a(lk)
843  lk = lk + 1
844  endif
845  if (j.lt.ny) then
846  stencil(5) = a(lk)
847  lk = lk + 1
848  end if
849  if (k.lt.nz) stencil(7) = a(lk)
850 c
851 c first the ia pointer -- points to the beginning of each row
852  ia(node) = iedge
853 c
854 c move the values from the local stencil to the new matrix
855 c
856 c the neighbor on the left (west)
857  if (i.gt.lx) then
858  ja(iedge)=node-kx
859  a(iedge) =stencil(2)
860  iedge=iedge + 1
861  end if
862 c the neighbor below (south)
863  if (j.gt.ly) then
864  ja(iedge)=node-ky
865  a(iedge)=stencil(4)
866  iedge=iedge + 1
867  end if
868 c the neighbor in the front
869  if (k.gt.lz) then
870  ja(iedge)=node-kz
871  a(iedge)=stencil(6)
872  iedge=iedge + 1
873  endif
874 c center node (itself)
875  ja(iedge) = node
876  iau(node) = iedge
877  a(iedge) = stencil(1)
878  iedge = iedge + 1
879 c the neighbor to the right (east)
880  if (i.lt.ux) then
881  ja(iedge)=node+kx
882  a(iedge)=stencil(3)
883  iedge=iedge + 1
884  end if
885 c the neighbor above (north)
886  if (j.lt.uy) then
887  ja(iedge)=node+ky
888  a(iedge)=stencil(5)
889  iedge=iedge + 1
890  end if
891 c the neighbor at the back
892  if (k.lt.uz) then
893  ja(iedge)=node+kz
894  a(iedge)=stencil(7)
895  iedge=iedge + 1
896  end if
897 c the right-hand side
898  rhs(node) = rhs(nbnode)
899 c------next node -------------------------
900  node=node+1
901 c
902  82 continue
903  81 continue
904  80 continue
905 c
906  ia(node) = iedge
907 c
908 c the number of nodes in the final matrix is stored in n
909 c
910  n = node - 1
911  return
912 c-----------------------------------------------------------------------
913  end
914 c-----end of fdreduce-----------------------------------------------------
915 c-----------------------------------------------------------------------
916  subroutine fdaddbc(nx,ny,nz,a,ja,ia,iau,rhs,al,h)
917  integer nx, ny, nz, ia(nx*ny*nz), ja(7*nx*ny*nz), iau(nx*ny*nz)
918  real*8 h, al(6), a(7*nx*ny*nz), rhs(nx*ny*nz)
919 c-----------------------------------------------------------------------
920 c This subroutine will add the boundary condition to the linear system
921 c consutructed without considering the boundary conditions
922 c
923 c The Boundary condition is specified in the following form:
924 c du
925 c alpha -- + beta u = gamma
926 c dn
927 c Alpha is stored in array AL. The six side of the boundary appares
928 c in AL in the following order: left(west), right(east), bottom(south),
929 c top(north), front, back(rear). (see also the illustration in gen57pt)
930 c Beta and gamma appears as the functions, betfun and gamfun.
931 c They have the following prototype
932 c
933 c real*8 function xxxfun(x, y, z)
934 c real*8 x, y, z
935 c
936 c where x, y, z are vales in the range of [0, 1][0, (ny-1)*h]
937 c [0, (nz-1)*h]
938 c
939 c At the corners or boundary lines, the boundary conditions are applied
940 c in the follow order:
941 c 1) if one side is Dirichlet boundary condition, the Dirichlet boundary
942 c condition is used;
943 c 2) if more than one sides are Dirichlet, the Direichlet condition
944 c specified for X direction boundary will overwrite the one specified
945 c for Y direction boundary which in turn has priority over Z
946 c direction boundaries.
947 c 3) when all sides are non-Dirichlet, the average values are used.
948 c-----------------------------------------------------------------------
949 c some constants
950 c
951  real*8 half,zero,one,two
952  parameter(half=0.5d0,zero=0.0d0,one=1.0d0,two=2.0d0)
953 c
954 c local variables
955 c
956  character*2 side
957  integer i,j,k,kx,ky,kz,node,nbr,ly,uy,lx,ux
958  real*8 coeff, ctr, hhalf, x, y, z
959  real*8 afun, bfun, cfun, dfun, efun, ffun, gfun, hfun
960  external afun, bfun, cfun, dfun, efun, ffun, gfun, hfun
961  real*8 betfun, gamfun
962  integer lctcsr
963  external lctcsr, betfun, gamfun
964 c
965  hhalf = half * h
966  kx = 1
967  ky = nx
968  kz = nx*ny
969 c
970 c In 3-D case, we need to go through all 6 faces one by one. If
971 c the actual dimension is lower, test on ny is performed first.
972 c If ny is less or equals to 1, then the value of nz is not
973 c checked.
974 c-----
975 c the surface on the left (west) side
976 c Concentrate on the contribution from the derivatives related to x,
977 c The terms with derivative of x was assumed to be:
978 c
979 c a(3/2,j,k)*[u(2,j,k)-u(1,j,k)] + a(1/2,j,k)*[u(0,j,k)-u(1,j,k)] +
980 c h*d(1,j,k)*[u(2,j,k)-u(0,j,k)]/2
981 c
982 c But they actually are:
983 c
984 c 2*{a(3/2,j,k)*[u(2,j,k)-u(1,j,k)] -
985 c h*a(1,j,k)*[beta*u(1,j,k)-gamma]/alpha]} +
986 c h*h*d(1,j,k)*[beta*u(1,j,k)-gamma]/alpha
987 c
988 c Therefore, in terms of local stencil the right neighbor of a node
989 c should be changed to 2*a(3/2,j,k),
990 c The matrix never contains the left neighbor on this border, nothing
991 c needs to be done about it.
992 c The following terms should be added to the center stencil:
993 c -a(3/2,j,k) + a(1/2,j,k) + [h*d(1,j,k)-2*a(1,j,k)]*h*beta/alpha
994 c
995 c And these terms should be added to the corresponding right-hand side
996 c [h*d(1,j,k)-2*a(1,j,k)]*h*gamma/alpha
997 c
998 c Obviously, the formula do not apply for the Dirichlet Boundary
999 c Condition, where alpha will be zero. In that case, we simply set
1000 c all the elements in the corresponding row to zero(0), then let
1001 c the diagonal element be beta, and the right-hand side be gamma.
1002 c Thus the value of u at that point will be set. Later on point
1003 c like this will be removed from the matrix, since they are of
1004 c know value before solving the system.(not done in this subroutine)
1005 c
1006  x = zero
1007  side = 'x1'
1008  do 20 k = 1, nz
1009  z = (k-1)*h
1010  do 21 j = 1, ny
1011  y = (j-1)*h
1012  node = 1+(j-1)*ky+(k-1)*kz
1013 c
1014 c check to see if it's Dirichlet Boundary condition here
1015 c
1016  if (al(1) .eq. zero) then
1017  call clrow(node, a, ja, ia)
1018  a(iau(node)) = betfun(side,x,y,z)
1019  rhs(node) = gamfun(side,x,y,z)
1020  else
1021 c
1022 c compute the terms formulated above to modify the matrix.
1023 c
1024 c the right neighbor is stroed in nbr'th posiiton in the a
1025  nbr = lctcsr(node, node+kx, ja, ia)
1026 c
1027  coeff = two*afun(x,y,z)
1028  ctr = (h*dfun(x,y,z) - coeff)*h/al(1)
1029  rhs(node) = rhs(node) + ctr * gamfun(side,x,y,z)
1030  ctr = afun(x-hhalf,y,z) + ctr * betfun(side,x,y,z)
1031  coeff = afun(x+hhalf,y,z)
1032  a(iau(node)) = a(iau(node)) - coeff + ctr
1033  a(nbr) = two*coeff
1034  end if
1035  21 continue
1036  20 continue
1037 c
1038 c the right (east) side boudary, similarly, the contirbution from
1039 c the terms containing the derivatives of x were assumed to be
1040 c
1041 c a(nx+1/2,j,k)*[u(nx+1,j,k)-u(nx,j,k)] +
1042 c a(nx-1/2,j,k)*[u(nx-1,j,k)-u(nx,j,k)] +
1043 c d(nx,j,k)*[u(nx+1,j,k)-u(nx-1,j,k)]*h/2
1044 c
1045 c Actualy they are:
1046 c
1047 c 2*{h*a(nx,j,k)*[gamma-beta*u(nx,j,k)]/alpha +
1048 c a(nx-1/2,j,k)*[u(nx-1,j,k)-u(nx,j,k)]} +
1049 c h*h*d(nx,j,k)*[gamma-beta*u(nx,j,k)]/alpha
1050 c
1051 c The left stencil has to be set to 2*a(nx-1/2,j,k)
1052 c
1053 c The following terms have to be added to the center stencil:
1054 c
1055 c -a(nx-1/2,j,k)+a(nx+1/2,j,k)-[2*a(nx,j,k)+h*d(nx,j,k)]*beta/alpha
1056 c
1057 c The following terms have to be added to the right-hand side:
1058 c
1059 c -[2*a(nx,j,k)+h*d(nx,j,k)]*h*gamma/alpha
1060 c
1061  x = one
1062  side = 'x2'
1063  do 22 k = 1, nz
1064  z = (k-1)*h
1065  do 23 j = 1, ny
1066  y = (j-1)*h
1067  node = (k-1)*kz + j*ky
1068 c
1069  if (al(2) .eq. zero) then
1070  call clrow(node, a, ja, ia)
1071  a(iau(node)) = betfun(side,x,y,z)
1072  rhs(node) = gamfun(side,x,y,z)
1073  else
1074  nbr = lctcsr(node, node-kx, ja, ia)
1075 c
1076  coeff = two*afun(x,y,z)
1077  ctr = (coeff + h*dfun(x,y,z))*h/al(2)
1078  rhs(node) = rhs(node) - ctr * gamfun(side,x,y,z)
1079  ctr = afun(x+hhalf,y,z) - ctr * betfun(side,x,y,z)
1080  coeff = afun(x-hhalf,y,z)
1081  a(iau(node)) = a(iau(node)) - coeff + ctr
1082  a(nbr) = two*coeff
1083  end if
1084  23 continue
1085  22 continue
1086 c
1087 c If only one dimension, return now
1088 c
1089  if (ny .le. 1) return
1090 c
1091 c the bottom (south) side suface, This similar to the situation
1092 c with the left side, except all the function and realted variation
1093 c should be on the y.
1094 c
1095 c These two block if statment here is to resolve the possible conflict
1096 c of assign the boundary value differently by different side of the
1097 c Dirichlet Boundary Conditions. They ensure that the edges that have
1098 c be assigned a specific value will not be reassigned.
1099 c
1100  if (al(1) .eq. zero) then
1101  lx = 2
1102  else
1103  lx = 1
1104  end if
1105  if (al(2) .eq. zero) then
1106  ux = nx-1
1107  else
1108  ux = nx
1109  end if
1110  y = zero
1111  side = 'y1'
1112  do 24 k = 1, nz
1113  z = (k-1)*h
1114  do 25 i = lx, ux
1115  x = (i-1)*h
1116  node = i + (k-1)*kz
1117 c
1118  if (al(3) .eq. zero) then
1119  call clrow(node, a, ja, ia)
1120  a(iau(node)) = betfun(side,x,y,z)
1121  rhs(node) = gamfun(side,x,y,z)
1122  else
1123  nbr = lctcsr(node, node+ky, ja, ia)
1124 c
1125  coeff = two*bfun(x,y,z)
1126  ctr = (h*efun(x,y,z) - coeff)*h/al(3)
1127  rhs(node) = rhs(node) + ctr * gamfun(side,x,y,z)
1128  ctr = bfun(x,y-hhalf,z) + ctr * betfun(side,x,y,z)
1129  coeff = bfun(x,y+hhalf,z)
1130  a(iau(node)) = a(iau(node)) - coeff + ctr
1131  a(nbr) = two*coeff
1132  end if
1133  25 continue
1134  24 continue
1135 c
1136 c The top (north) side, similar to the right side
1137 c
1138  y = (ny-1) * h
1139  side = 'y2'
1140  do 26 k = 1, nz
1141  z = (k-1)*h
1142  do 27 i = lx, ux
1143  x = (i-1)*h
1144  node = (k-1)*kz+(ny-1)*ky + i
1145 c
1146  if (al(4) .eq. zero) then
1147  call clrow(node, a, ja, ia)
1148  a(iau(node)) = betfun(side,x,y,z)
1149  rhs(node) = gamfun(side,x,y,z)
1150  else
1151  nbr = lctcsr(node, node-ky, ja, ia)
1152 c
1153  coeff = two*bfun(x,y,z)
1154  ctr = (coeff + h*efun(x,y,z))*h/al(4)
1155  rhs(node) = rhs(node) - ctr * gamfun(side,x,y,z)
1156  ctr = bfun(x,y+hhalf,z) - ctr * betfun(side,x,y,z)
1157  coeff = bfun(x,y-hhalf,z)
1158  a(iau(node)) = a(iau(node)) - coeff + ctr
1159  a(nbr) = two*coeff
1160  end if
1161  27 continue
1162  26 continue
1163 c
1164 c If only has two dimesion to work on, return now
1165 c
1166  if (nz .le. 1) return
1167 c
1168 c The front side boundary
1169 c
1170 c If the edges of the surface has been decided by Dirichlet Boundary
1171 c Condition, then leave them alone.
1172 c
1173  if (al(3) .eq. zero) then
1174  ly = 2
1175  else
1176  ly = 1
1177  end if
1178  if (al(4) .eq. zero) then
1179  uy = ny-1
1180  else
1181  uy = ny
1182  end if
1183 c
1184  z = zero
1185  side = 'z1'
1186  do 28 j = ly, uy
1187  y = (j-1)*h
1188  do 29 i = lx, ux
1189  x = (i-1)*h
1190  node = i + (j-1)*ky
1191 c
1192  if (al(5) .eq. zero) then
1193  call clrow(node, a, ja, ia)
1194  a(iau(node)) = betfun(side,x,y,z)
1195  rhs(node) = gamfun(side,x,y,z)
1196  else
1197  nbr = lctcsr(node, node+kz, ja, ia)
1198 c
1199  coeff = two*cfun(x,y,z)
1200  ctr = (h*ffun(x,y,z) - coeff)*h/al(5)
1201  rhs(node) = rhs(node) + ctr * gamfun(side,x,y,z)
1202  ctr = cfun(x,y,z-hhalf) + ctr * betfun(side,x,y,z)
1203  coeff = cfun(x,y,z+hhalf)
1204  a(iau(node)) = a(iau(node)) - coeff + ctr
1205  a(nbr) = two*coeff
1206  end if
1207  29 continue
1208  28 continue
1209 c
1210 c Similiarly for the top side of the boundary suface
1211 c
1212  z = (nz - 1) * h
1213  side = 'z2'
1214  do 30 j = ly, uy
1215  y = (j-1)*h
1216  do 31 i = lx, ux
1217  x = (i-1)*h
1218  node = (nz-1)*kz + (j-1)*ky + i
1219 c
1220  if (al(6) .eq. zero) then
1221  call clrow(node, a, ja, ia)
1222  a(iau(node)) = betfun(side,x,y,z)
1223  rhs(node) = gamfun(side,x,y,z)
1224  else
1225  nbr = lctcsr(node, node-kz, ja, ia)
1226 c
1227  coeff = two*cfun(x,y,z)
1228  ctr = (coeff + h*ffun(x,y,z))*h/al(6)
1229  rhs(node) = rhs(node) - ctr * gamfun(side,x,y,z)
1230  ctr = cfun(x,y,z+hhalf) - ctr * betfun(side,x,y,z)
1231  coeff = cfun(x,y,z-hhalf)
1232  a(iau(node)) = a(iau(node)) - coeff + ctr
1233  a(nbr) = two*coeff
1234  end if
1235  31 continue
1236  30 continue
1237 c
1238 c all set
1239 c
1240  return
1241 c-----------------------------------------------------------------------
1242  end
1243 c-----end of fdaddbc----------------------------------------------------
1244 c-----------------------------------------------------------------------
1245  subroutine clrow(i, a, ja, ia)
1246  integer i, ja(*), ia(*), k
1247  real *8 a(*)
1248 c-----------------------------------------------------------------------
1249 c clear the row i to all zero, but still keep the structure of the
1250 c CSR matrix
1251 c-----------------------------------------------------------------------
1252  do 10 k = ia(i), ia(i+1)-1
1253  a(k) = 0.0d0
1254  10 continue
1255 c
1256  return
1257 c-----end of clrow------------------------------------------------------
1258  end
1259 c-----------------------------------------------------------------------
1260  function lctcsr(i,j,ja,ia)
1261  integer lctcsr, i, j, ja(*), ia(*), k
1262 c-----------------------------------------------------------------------
1263 c locate the position of a matrix element in a CSR format
1264 c returns -1 if the desired element is zero
1265 c-----------------------------------------------------------------------
1266  lctcsr = -1
1267  k = ia(i)
1268  10 if (k .lt. ia(i+1) .and. (lctcsr .eq. -1)) then
1269  if (ja(k) .eq. j) lctcsr = k
1270  k = k + 1
1271  goto 10
1272  end if
1273 c
1274  return
1275 c-----------------------------------------------------------------------
1276  end
1277 c-----end of lctcsr-----------------------------------------------------
1278 
1279  subroutine csrcsc (n,job,ipos,a,ja,ia,ao,jao,iao)
1280  integer ia(n+1),iao(n+1),ja(*),jao(*)
1281  real*8 a(*),ao(*)
1282 c-----------------------------------------------------------------------
1283 c Compressed Sparse Row to Compressed Sparse Column
1284 c
1285 c (transposition operation) Not in place.
1286 c-----------------------------------------------------------------------
1287 c -- not in place --
1288 c this subroutine transposes a matrix stored in a, ja, ia format.
1289 c ---------------
1290 c on entry:
1291 c----------
1292 c n = dimension of A.
1293 c job = integer to indicate whether to fill the values (job.eq.1) of the
1294 c matrix ao or only the pattern., i.e.,ia, and ja (job .ne.1)
1295 c
1296 c ipos = starting position in ao, jao of the transposed matrix.
1297 c the iao array takes this into account (thus iao(1) is set to ipos.)
1298 c Note: this may be useful if one needs to append the data structure
1299 c of the transpose to that of A. In this case use for example
1300 c call csrcsc (n,1,ia(n+1),a,ja,ia,a,ja,ia(n+2))
1301 c for any other normal usage, enter ipos=1.
1302 c a = real array of length nnz (nnz=number of nonzero elements in input
1303 c matrix) containing the nonzero elements.
1304 c ja = integer array of length nnz containing the column positions
1305 c of the corresponding elements in a.
1306 c ia = integer of size n+1. ia(k) contains the position in a, ja of
1307 c the beginning of the k-th row.
1308 c
1309 c on return:
1310 c ----------
1311 c output arguments:
1312 c ao = real array of size nzz containing the "a" part of the transpose
1313 c jao = integer array of size nnz containing the column indices.
1314 c iao = integer array of size n+1 containing the "ia" index array of
1315 c the transpose.
1316 c
1317 c-----------------------------------------------------------------------
1318  call csrcsc2(n,n,job,ipos,a,ja,ia,ao,jao,iao)
1319  end
1320 c-----------------------------------------------------------------------
1321  subroutine csrcsc2 (n,n2,job,ipos,a,ja,ia,ao,jao,iao)
1322  integer ia(n+1),iao(n2+1),ja(*),jao(*)
1323  real*8 a(*),ao(*)
1324 c-----------------------------------------------------------------------
1325 c Compressed Sparse Row to Compressed Sparse Column
1326 c
1327 c (transposition operation) Not in place.
1328 c-----------------------------------------------------------------------
1329 c Rectangular version. n is number of rows of CSR matrix,
1330 c n2 (input) is number of columns of CSC matrix.
1331 c-----------------------------------------------------------------------
1332 c -- not in place --
1333 c this subroutine transposes a matrix stored in a, ja, ia format.
1334 c ---------------
1335 c on entry:
1336 c----------
1337 c n = number of rows of CSR matrix.
1338 c n2 = number of columns of CSC matrix.
1339 c job = integer to indicate whether to fill the values (job.eq.1) of the
1340 c matrix ao or only the pattern., i.e.,ia, and ja (job .ne.1)
1341 c
1342 c ipos = starting position in ao, jao of the transposed matrix.
1343 c the iao array takes this into account (thus iao(1) is set to ipos.)
1344 c Note: this may be useful if one needs to append the data structure
1345 c of the transpose to that of A. In this case use for example
1346 c call csrcsc2 (n,n,1,ia(n+1),a,ja,ia,a,ja,ia(n+2))
1347 c for any other normal usage, enter ipos=1.
1348 c a = real array of length nnz (nnz=number of nonzero elements in input
1349 c matrix) containing the nonzero elements.
1350 c ja = integer array of length nnz containing the column positions
1351 c of the corresponding elements in a.
1352 c ia = integer of size n+1. ia(k) contains the position in a, ja of
1353 c the beginning of the k-th row.
1354 c
1355 c on return:
1356 c ----------
1357 c output arguments:
1358 c ao = real array of size nzz containing the "a" part of the transpose
1359 c jao = integer array of size nnz containing the column indices.
1360 c iao = integer array of size n+1 containing the "ia" index array of
1361 c the transpose.
1362 c
1363 c-----------------------------------------------------------------------
1364 c----------------- compute lengths of rows of transp(A) ----------------
1365  do 1 i=1,n2+1
1366  iao(i) = 0
1367  1 continue
1368  do 3 i=1, n
1369  do 2 k=ia(i), ia(i+1)-1
1370  j = ja(k)+1
1371  iao(j) = iao(j)+1
1372  2 continue
1373  3 continue
1374 c---------- compute pointers from lengths ------------------------------
1375  iao(1) = ipos
1376  do 4 i=1,n2
1377  iao(i+1) = iao(i) + iao(i+1)
1378  4 continue
1379 c--------------- now do the actual copying -----------------------------
1380  do 6 i=1,n
1381  do 62 k=ia(i),ia(i+1)-1
1382  j = ja(k)
1383  next = iao(j)
1384  if (job .eq. 1) ao(next) = a(k)
1385  jao(next) = i
1386  iao(j) = next+1
1387  62 continue
1388  6 continue
1389 c-------------------------- reshift iao and leave ----------------------
1390  do 7 i=n2,1,-1
1391  iao(i+1) = iao(i)
1392  7 continue
1393  iao(1) = ipos
1394 c--------------- end of csrcsc2 ----------------------------------------
1395 c-----------------------------------------------------------------------
1396  end
1397 
subroutine gen57pt(nx, ny, nz, al, mode, n, a, ja, ia, iau, rhs)
Definition: genmat.f:20
integer function lctcsr(i, j, ja, ia)
Definition: genmat.f:1260
void csrcsc(int OUTINDEX, const int nrow, const int ncol, int job, double *a, int *ja, int *ia, double *ao, int *jao, int *iao)
convert csr to csc Assume input csr is 0-based index output csc 0/1 index specified by OUTINDEX * ...
Definition: spmat.c:17
subroutine gen57bl(nx, ny, nz, nfree, na, n, a, ja, ia, iau, stencil)
Definition: genmat.f:351
real *8 function efun(x, y, z)
Definition: functns.f90:32
real *8 function betfun(side, x, y, z)
Definition: functns.f90:58
subroutine bfunbl(nfree, x, y, z, coeff)
Definition: functns.f90:94
subroutine csrcsc2(n, n2, job, ipos, a, ja, ia, ao, jao, iao)
Definition: genmat.f:1321
subroutine dfunbl(nfree, x, y, z, coeff)
Definition: functns.f90:116
subroutine ffunbl(nfree, x, y, z, coeff)
Definition: functns.f90:136
real *8 function ffun(x, y, z)
Definition: functns.f90:40
subroutine bsten(nx, ny, nz, kx, ky, kz, nfree, stencil, h)
Definition: genmat.f:536
subroutine clrow(i, a, ja, ia)
Definition: genmat.f:1245
subroutine getsten(nx, ny, nz, mode, kx, ky, kz, stencil, h, rhs)
Definition: genmat.f:248
real *8 function cfun(x, y, z)
Definition: functns.f90:18
real *8 function gamfun(side, x, y, z)
Definition: functns.f90:65
subroutine fdaddbc(nx, ny, nz, a, ja, ia, iau, rhs, al, h)
Definition: genmat.f:916
real *8 function gfun(x, y, z)
Definition: functns.f90:46
subroutine cfunbl(nfree, x, y, z, coeff)
Definition: functns.f90:105
real *8 function dfun(x, y, z)
Definition: functns.f90:24
real *8 function hfun(x, y, z)
Definition: functns.f90:52
subroutine gfunbl(nfree, x, y, z, coeff)
Definition: functns.f90:146
subroutine efunbl(nfree, x, y, z, coeff)
Definition: functns.f90:126
real *8 function afun(x, y, z)
Definition: functns.f90:6
subroutine afunbl(nfree, x, y, z, coeff)
Definition: functns.f90:83
real *8 function bfun(x, y, z)
Definition: functns.f90:12
subroutine fdreduce(nx, ny, nz, alpha, n, a, ja, ia, iau, rhs, stencil)
Definition: genmat.f:656