Line data Source code
1 0 : subroutine nrk(x,y,zk,ap,aq,rhs,bc,ii,kk,ka,kb,ki,nn,id,ucy,ea,
2 0 : .det,v)
3 : c
4 : c ***********************************************************************
5 : c
6 : c package for solving boundary value problems for ordinary differential
7 : c equations with optional eigenvalues. accepts one or two point boundary
8 : c conditions, and integral constraints.
9 : c
10 : c method: replaces differential equations by centred, second order
11 : c accuracy difference equations, and solves these by
12 : c newton-raphson-kantorovich iteration, using forward
13 : c elimination and back-substitution. see baker, moore and
14 : c spiegel, q. j. mech. appl. math., vol. 24, p. 391
15 : c (1971).
16 : c
17 : c package written by d. o. gough, institute of astronomy and
18 : c department of applied mathematics and theoretical phyisics,
19 : c university of cambridge, england.
20 : c
21 : c **********************************************************************
22 : c
23 : c newton-raphson-kantorovich setting up routine
24 : c
25 : c space required by common/work/ is
26 : c
27 : c (2*((ii+kk+1)*(6*ii+2*kk+1)+ii*(ii-ka)+ki*(ii+kk+2))
28 : c + (ii+kk+1-ka)*ii*nn)*4 bytes
29 : c
30 : c the programmer has the option of defining common/nrkchk/nrchk
31 : c and setting it to the number of single precision words that has been
32 : c reserved in labelled common/work/
33 : c if nrkwsp.gt.0, nrk checks that it is great enough to accommodate th
34 : c problem.
35 : c if the work space is too small or if nrkwsp.le.0, appropriate
36 : c diagnostics are written
37 : c
38 : c if an error is detected v(1) is set to zero before return
39 : c
40 : c
41 : c Double precision version.
42 : c ++++++++++++++++++++++++
43 : c
44 : c Dated 16/4/90.
45 : c
46 : c Note: this double precision version of the routine has same name
47 : c as single precision version, but is distinguished by its file name
48 : c
49 : implicit double precision (a-h,o-z)
50 : integer v
51 : dimension x(nn),y(id,nn),zk(kk),ea(id,3),ap(1),aq(1),v(ii)
52 : common/nrkchk/nrkwsp
53 : common/work/ w(2000)
54 : c
55 : c common defining standard input and output
56 : c
57 : common/cstdio/ istdin, istdou, istdpr, istder
58 : external rhs,bc
59 : data ichk/0/
60 : c
61 : c
62 : c line printer dsrn
63 : c
64 : c set dimension parameters
65 : c
66 0 : 20 ik=ii+kk
67 0 : ik1=ik+1
68 0 : ikka1=ik-ka+1
69 0 : ip=ii*ikka1
70 0 : ipn=ip*nn
71 0 : iiik=ii*ik
72 0 : iiik1=ii+ik1
73 0 : kid=ki
74 0 : if(ki.le.0)kid=1
75 : c
76 : c location of variables
77 : c note: a,d,p must be stored consecutively and in order
78 : c
79 0 : 30 lq=1
80 0 : lf=lq+ip
81 0 : lh=lf+ii+ii
82 0 : lfd=lh+kid
83 0 : lhd=lfd+iiik+iiik
84 0 : lg=lhd+ki*ik1
85 0 : lgp=lg+ik1
86 0 : la=lgp+2*ik1*ik1
87 0 : ld=la+ii*iiik1
88 0 : lp=ld+iiik+ii
89 : c
90 : c size of common work/
91 : c
92 0 : iwork=(lp-lq)+ ipn
93 : c
94 : c if nrkwsp has been set, check that enough space has been reserved
95 : c in /work/
96 : c
97 0 : 40 if(nrkwsp.gt.0)go to 42
98 0 : 41 if(ichk.ne.iwork.and.istdpr.gt.0) write(istdpr,1001)iwork
99 0 : ichk=iwork
100 0 : go to 50
101 0 : 42 if(nrkwsp.ge.iwork)go to 50
102 0 : 43 write(istdou,1002)nrkwsp,iwork
103 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdpr,1002)
104 0 : * nrkwsp,iwork
105 0 : v(1)=0
106 0 : return
107 : c
108 : c call newton-raphson-kantorovich routine
109 : c note:logically gd and gp are equivalent
110 : c g,gst,eb and gs,p share store to save space, consequently ik1*8 byte
111 : c have been reserved for g in /work/ to accommodate gst
112 : c
113 : 50 call nrke(x,y,zk,ap,aq,rhs,bc,ii,kk,ka,kb,ki,nn,id,ucy,ea,det,v
114 : .,w(lf),w(lfd),w(lg),w(lg),w(lgp),w(lgp),w(lp),w(la),w(ld),w(lh)
115 0 : .,w(lhd),w(lq),w(lp),w(lg),ik,ik1,ikka1,iiik1,ip,ipn,iw,kid)
116 0 : return
117 : c
118 : c diagnostic messages
119 : c
120 : 1001 format(1x,5('%'),3x,'size of common/work/ not checked by nrk',
121 : .i9,'*4bytes required ... execution continuing')
122 : 1002 format(1x,5('%'),i9,'*4 bytes in common/work/ insufficient for
123 : . nrk',i9,'*4 bytes required ... nrk execution terminated')
124 : c
125 : end
126 0 : subroutine nrke(x,y,zk,ap,aq,rhs,bc,ii,kk,ka,kb,ki,nn,id,ucy,ea,
127 0 : . det,v,f,fd,g,gst,gd,gp,gs,a,d,h,hd,q,p,eb,ik,
128 : . ik1,ikka1,iiik1,ip,ipn,iw,kid)
129 : c
130 : c
131 : c newton-raphson-kantorovich programme
132 : c ************************************
133 : c
134 : c if an error is detected, v(1) is set to zero before return
135 : c
136 : c
137 : implicit double precision (a-h,o-z)
138 : integer v
139 : dimension f(ii,2),fd(ii,ik,2),g(ik),gst(ik1),gd(ik1,ik1),
140 : . gp(ik1,ik1,2),gs(ik,ik1,2),a(ii,iiik1),d(ii,ik1),
141 : . h(kid),hd(kid,ik1),q(ii,ikka1),p(ipn)
142 : c
143 : dimension x(nn),y(id,nn),zk(kk),ea(id,3),
144 : . eb(ii),v(ii),ap(1),aq(1)
145 : c
146 : c
147 : external rhs,bc
148 : c note: rhs sets derivatives of h into d and not hd. therefore
149 : c it must assume that the first dimension of hd is the
150 : c same as that of fd
151 : c
152 : c
153 : c
154 : c
155 : c line printer dsrn
156 : c
157 : c set counting limits
158 0 : kab=ka+kb
159 0 : kap=ik-ki
160 0 : ika=ii-ka
161 0 : ikka=ik-ka
162 0 : kp=kap-kab
163 0 : kip=ki+kp
164 : c
165 0 : 3 n1=nn-1
166 0 : i1=ii+1
167 0 : k1=kk+1
168 0 : kab1=kab+1
169 0 : kap1=kap+1
170 0 : ka1=ka+1
171 0 : ika1=ika+1
172 0 : kaik=ka*ik1
173 0 : ikka1k=ikka1*ik1
174 0 : ii2=ii*ii
175 0 : ikka1i=ikka1*ii
176 0 : ikkaik=ikka*ik1
177 : c
178 : c
179 : c compatibility test
180 0 : if(ka.le.ii.and.kp.ge.0.and.ka.ge.0.and.kb.ge.0.and.ki.ge.0.and.
181 0 : . kk.ge.0) go to 10
182 0 : write(istdou,1000) ii,kk,ka,kb,ki
183 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdpr,1000)
184 0 : * ii,kk,ka,kb,ki
185 0 : v(1)=0
186 0 : return
187 : 10 continue
188 : c
189 : c
190 : c set v if necessary
191 0 : 20 iv=0
192 0 : do 21 i=1,ii
193 0 : if(v(i).gt.0) go to 21
194 0 : iv=1
195 0 : v(i)=i
196 0 : if(istdpr.gt.0) write(istdpr,1100) i,i
197 0 : 21 continue
198 0 : if(iv.eq.1.and.istdpr.gt.0) write(istdpr,1101) (v(i), i=1,ii)
199 : c
200 : c set range of independent variable
201 0 : 25 r=x(nn)-x(1)
202 0 : if(r.eq.0.) go to 407
203 : c
204 : c
205 : c **********
206 : c conditions solely at first boundary
207 : c
208 : c empty derivative matrices
209 0 : 30 do 34 j=1,ik
210 0 : do 31 i=1,ii
211 0 : 31 fd(i,j,1)=0.0
212 0 : do 32 k=1,2
213 0 : do 32 i=1,kap
214 0 : 32 gs(i,j,k)=0.0
215 0 : if(ki.eq.0) go to 34
216 0 : do 33 ig=1,ki
217 0 : 33 d(ig,j)=0.0
218 0 : 34 continue
219 : c
220 : c
221 : c boundary conditions and equations at first boundary
222 0 : 35 ia=ik
223 0 : ib=ik1
224 0 : in=ii
225 0 : i=1
226 0 : 36 call bc(x(1),x(nn),y,y(1,nn),zk,ap,aq,g,gs,ia,ib,nn)
227 0 : call rhs(x(1),y,zk,ap,aq,f,fd,h,d,in,i)
228 : c note that although the dimensions of d might not be adequate to
229 : c store hd, in that event there is always sufficient unused space
230 : c at the beginning of p to accommodate the overflow
231 : c
232 : c set boundary derivative matrix
233 0 : 40 do 43 i=1,ik
234 0 : iv=i
235 0 : if(i.le.ii) iv=v(i)
236 0 : do 41 j=1,kap
237 0 : do 41 k=1,2
238 0 : 41 gp(j,i,k)=gs(j,iv,k)
239 : c (ensure that first and second b.c. matrices are loaded into gd)
240 0 : if(kab.eq.0) go to 43
241 0 : do 42 j=1,kab
242 0 : 42 gd(j,i)=gd(j,i)+gp(j,i,2)
243 0 : 43 gd(i,ik1)=g(i)
244 : c (ensure that eigenvalue dependence is loaded into gd)
245 0 : if(kp.eq.0.or.kk.eq.0) go to 45
246 0 : do 44 is=kab1,kap
247 0 : do 44 j=i1,ik
248 0 : 44 gd(is,j)=gd(is,j)+gp(is,j,2)
249 : c
250 : c first inversion (for equation 5)
251 0 : 45 detsgn=1.
252 0 : if(ka.ne.0) go to 46
253 0 : det=0.
254 0 : go to 48
255 0 : 46 call leqsd(gd,gd(1,ka1),ka,ikka1,ik1,ik1,err,kaik,ikka1k)
256 0 : if(err) 47,400,48
257 0 : 47 detsgn=-detsgn
258 0 : 48 det=log10( abs(err))
259 : c note that now gd(ial,nu)=-h(ial,nu,1) in equation 5
260 : c
261 : c initial integration coefficients
262 0 : hn=(x(2)-x(1))/6.e0
263 0 : hn1=(x(3)-x(2))/6.e0
264 0 : if(ki.eq.0) go to 50
265 0 : wa=hn+hn1
266 0 : wb=hn-hn1
267 0 : wd=wa/hn
268 0 : ah=wd*(hn+wb)
269 0 : ah1=wa*wa*wd/hn1
270 0 : ah2=wa*(hn1-wb)/hn1
271 : c set integral constraint matrix
272 0 : do 49 ig=1,ki
273 0 : hd(ig,ik1)=h(ig)
274 0 : do 49 i=1,ik
275 0 : iv=i
276 0 : if(i.le.ii) iv=v(i)
277 0 : 49 hd(ig,i)=d(ig,iv)
278 : c
279 : c
280 : c compute p and gamma (gd) at first mesh point
281 : c and store p in real*8 q for computation of e in loops 103 and 105
282 0 : 50 if(ka.eq.0) go to 55
283 0 : do 54 nu=ka1,ik1
284 0 : j=(nu-ka-1)*ii+ika
285 0 : do 51 ial=1,ka
286 0 : wd=-gd(ial,nu)
287 0 : q(ika+ial,nu-ka)=wd
288 0 : 51 p(ial+j)=wd
289 0 : if(kp.eq.0) go to 54
290 0 : do 53 is=kab1,kap
291 0 : wd=0.0
292 0 : do 52 ial=1,ka
293 0 : 52 wd=wd-gd(is,ial)*gd(ial,nu)
294 0 : 53 gd(is,nu)=wd+gd(is,nu)
295 0 : 54 continue
296 : c contribution at first boundary to integral constraints
297 0 : 55 if(ki.eq.0) go to 58
298 0 : do 57 nu=ka1,ik1
299 0 : do 57 ig=1,ki
300 0 : wd=0.0
301 0 : if(ka.eq.0) go to 57
302 0 : do 56 ial=1,ka
303 0 : 56 wd=wd-hd(ig,ial)*gd(ial,nu)
304 0 : 57 gd(ig+kap,nu)=ah*(wd+hd(ig,nu))
305 : 58 continue
306 : c
307 : c
308 : c
309 : c **********
310 : c beginning of preliminary outer loop
311 : c
312 : c set storage indices
313 0 : 70 in=1
314 0 : in1=2
315 : c
316 0 : 80 do 130 n=1,n1
317 0 : np1=n+1
318 0 : dx=3.e0*hn
319 : c
320 : c check that independent variable is strictly monotonic
321 0 : if(r*dx.le.0.) go to 408
322 : c
323 : c compute integration coefficients
324 0 : hn=hn1
325 0 : if(n+3.gt.nn) go to 81
326 0 : hn1=(x(n+3)-x(n+2))/6.e0
327 0 : 81 if(ki.eq.0) go to 90
328 0 : if(n.eq.n1) go to 85
329 0 : if(in.eq.1) go to 84
330 : c check whether np1=nn-1 when nn is even
331 0 : if(np1.eq.n1) go to 84
332 0 : wa=hn+hn1
333 0 : wb=hn-hn1
334 0 : wd=wa/hn
335 0 : ah=wd*(hn+wb)+ah2
336 0 : ah2=wa/hn1
337 0 : ah1=wa*wd*ah2
338 0 : ah2=ah2*(hn1-wb)
339 0 : go to 90
340 : c integration coefficients at np1 when np1 is even
341 : c and nn-1 when nn is even
342 0 : 84 ah=ah1
343 : c check whether np1=nn-2
344 0 : if(n.ne.nn-3) go to 90
345 : c integration coefficients at nn-2,nn-1 and nn when nn is even
346 0 : wa=hn+hn1
347 0 : wb=3.e0*hn+hn1
348 0 : ah=ah-hn1*hn1*hn1/(hn*wa)
349 0 : ah1=hn1*wb/hn+ah2
350 0 : ah2=hn1*(wb+hn1)/wa
351 0 : go to 90
352 : c integration coefficient at nn
353 0 : 85 ah=ah2
354 : c
355 : c equations at mesh point np1
356 0 : 90 do 92 i=1,ik
357 0 : do 91 j=1,ii
358 0 : 91 fd(j,i,in1)=0.e0
359 0 : do 92 ig=1,ki
360 0 : 92 d(ig,i)=0.e0
361 0 : ia=ii
362 0 : i=np1
363 0 : call rhs(x(np1),y(1,np1),zk,ap,aq,f(1,in1),fd(1,1,in1),h,d,ia,i)
364 0 : do 93 i=1,ii
365 0 : 93 gst(i)=y(v(i),n)-y(v(i),np1)
366 0 : 96 if(ki.eq.0) go to 100
367 0 : do 97 i=1,ik
368 0 : iv=i
369 0 : if(i.le.ii) iv=v(i)
370 0 : do 97 ig=1,ki
371 0 : 97 hd(ig,i)=d(ig,iv)
372 : c
373 : c beginning of loop to construct e matrix
374 0 : 100 do 109 i=1,ii
375 : c
376 : c construct a and b matrices (b is loaded into d)
377 0 : do 101 j=1,ii
378 0 : a(i,j)=dx*fd(v(i),v(j),in)
379 0 : 101 d(i,j)=dx*fd(v(i),v(j),in1)
380 0 : a(i,i)=a(i,i)+1.0e0
381 0 : d(i,i)=d(i,i)-1.0e0
382 0 : do 102 k= i1,ik
383 0 : 102 d(i,k)=dx*(fd(v(i),k,in)+fd(v(i),k,in1))
384 0 : d(i,ik1)=dx*(f(v(i),in)+f(v(i),in1))+gst(i)
385 : c construct d matrix (b block is already loaded)
386 0 : if(ka.eq.0) go to 107
387 0 : do 104 m=i1,ik1
388 0 : wd=0.0
389 0 : do 103 ial=1,ka
390 0 : 103 wd=wd+a(i,ial)*q(ika+ial,m-ka)
391 0 : 104 d(i,m)=d(i,m)+wd
392 : c
393 : c construct beginning of e matrix and set into end of a
394 0 : if(ika.eq.0) go to 107
395 0 : do 106 j=1,ika
396 0 : wd=0.0
397 0 : do 105 ial=1,ka
398 0 : 105 wd=wd+a(i,ial)*q(ika+ial,j)
399 0 : 106 a(i,j+ka)=wd+a(i,j+ka)
400 : 107 continue
401 : c
402 : c douglas is playing around here and we need to dort it out ]]
403 : c
404 0 : do 108 m=1,ik1
405 0 : 108 a(i,ii+m)=d(i,m)
406 : c
407 : c in 110 and 111,121 loops a(i,ii+m) will be used in place of d(i,m)
408 0 : 109 continue
409 : c end of loop to construct e matrix
410 : c
411 : c
412 : c compute p matrix
413 : c and store current value in real*8 q for computation of e
414 : c in loops 103 and 105
415 : c
416 0 : 110 call leqsd(a(1,ka1),d(1,ka1),ii,ikka1,ii,ii,err,ii2,ikka1i)
417 0 : if(err) 111,403,112
418 0 : 111 detsgn=-detsgn
419 0 : 112 k=n*ip
420 0 : det=det+log10( abs(err))
421 0 : do 115 nu=ka1,ik1
422 0 : j=(nu-ka-1)*ii+k
423 0 : do 115 i=1,ii
424 0 : wd=-d(i,nu)
425 0 : q(i,nu-ka)=wd
426 0 : 115 p(i+j)=wd
427 : c
428 : c compute gamma matrix
429 0 : 120 if(kip.eq.0.or.ika.eq.0) go to 125
430 0 : do 124 is=kab1,ik
431 0 : do 123 nu=ka1,ik1
432 0 : wd=0.0
433 0 : do 121 ib=ka1,ii
434 0 : 121 wd=wd-gd(is,ib)*d(ib-ka,nu)
435 0 : if(nu.gt.ii) go to 122
436 0 : gst(nu-ka)=wd
437 0 : go to 123
438 0 : 122 gst(nu-ka)=wd+gd(is,nu)
439 0 : 123 continue
440 0 : do 124 nu=ka1,ik1
441 0 : 124 gd(is,nu)=gst(nu-ka)
442 : c integral constraints at point np1
443 0 : 125 if(ki.eq.0) go to 129
444 0 : do 127 ig=1,ki
445 0 : hd(ig,ik1)=h(ig)
446 0 : do 127 nu=ka1,ik1
447 0 : wd=0.0
448 0 : if(ka.eq.0) go to 127
449 0 : do 126 ial=1,ka
450 0 : 126 wd=wd-hd(ig,ial)*d(ika+ial,nu)
451 0 : 127 gd(ig+kap,nu)=gd(ig+kap,nu)+ah*(wd+hd(ig,nu))
452 : c
453 : c reset storage indices
454 0 : 129 i=in
455 0 : in=in1
456 0 : in1=i
457 0 : 130 continue
458 : c
459 : c end of preliminary outer loop
460 : c **********
461 : c
462 : c
463 : c **********
464 : c remaining boundary conditions
465 : c
466 0 : 200 if(ka.eq.ik) go to 233
467 : c compute coefficients for equation 12a and load into end of gd
468 0 : if(kp.eq.0) go to 210
469 0 : do 205 is=kab1,kap
470 0 : do 202 nu=ka1,ik1
471 0 : wd=0.0
472 0 : do 201 ial=1,ka1
473 0 : 201 wd=wd-gp(is,ial,2)*d(ika+ial,nu)
474 0 : 202 gd(is,nu)=gd(is,nu)+wd
475 0 : do 203 ia=ka1,ii
476 0 : 203 gd(is,ia)=gd(is,ia)+gp(is,ia,2)
477 0 : 205 continue
478 : 210 continue
479 : c
480 : c
481 : c compute coefficients for equation 13 and load into middle of gd
482 0 : 220 if(ka.eq.0.or.kb.eq.0) go to 230
483 0 : do 223 m=ka1,kab
484 0 : do 222 nu=ka1,ik1
485 0 : wd=0.0
486 0 : do 221 ial=1,ka
487 0 : 221 wd=wd-gd(m,ial)*d(ika+ial,nu)
488 0 : 222 gd(m,nu)=wd+gd(m,nu)
489 0 : 223 continue
490 : c
491 : c solve equations 12a and 13 (answer has wrong sign)
492 0 : 230 if(ikka.eq.0) go to 233
493 : 231 call leqsd(gd(ka1,ka1),gd(ka1,ik1),ikka,1,ik1,ikka,err,ikkaik,
494 0 : . ikka)
495 0 : if(err) 233,404,234
496 0 : 233 detsgn=-detsgn
497 0 : 234 det=det+log10( abs(err))
498 0 : gd(ik1,ik1)=-1.0e0
499 0 : gd(ik1,1)=-1.0e0
500 : c (note that the first element of the second half of gp may now
501 : c have been overwritten)
502 : c
503 : c
504 : c
505 : c **********
506 : c iterated solution
507 : c
508 : c empty ea,eb
509 0 : 300 do 301 i=1,ii
510 0 : eb(i)=0.0
511 0 : do 301 j=1,3
512 0 : 301 ea(i,j)=0.0
513 : c
514 : c solution at second boundary
515 0 : 310 if(kk.eq.0) go to 313
516 0 : do 312 k=1,kk
517 0 : wd=gd(ii+k,ik1)
518 0 : gd(ii+k,1)=wd
519 0 : zk(k)=zk(k)-ucy*wd
520 0 : 312 continue
521 0 : 313 if(ka.ge.ii) go to 320
522 0 : do 318 ia=ka1,ii
523 0 : wa= abs(gd(ia,ik1))
524 0 : wd=y(v(ia),nn)
525 0 : 315 wd=wd-ucy*gd(ia,ik1)
526 0 : wb= abs(wd)
527 0 : ea(v(ia),1)=wa
528 0 : eb(ia)=wb
529 0 : if(wb.lt.1.0e-10) go to 316
530 0 : ea(v(ia),2)=wa/wb
531 0 : ea(v(ia),3)= float(nn)
532 : 316 continue
533 0 : y(v(ia),nn)=wd
534 0 : 318 continue
535 : c
536 : c remainder of solution (beginning of second outer loop)
537 : c
538 : c set storage indices
539 0 : 320 in=1
540 0 : in1=ik1
541 : c
542 0 : do 330 m=1,nn
543 0 : n=nn+1-m
544 0 : np1=ip*(n-1)
545 0 : do 328 i=1,ii
546 0 : if(i.gt.ika) go to 321
547 0 : k=n-1
548 0 : if(k.lt.1) go to 328
549 0 : j=i+ka
550 0 : go to 322
551 0 : 321 j=i-ika
552 0 : k=n
553 0 : 322 wd=0.0
554 0 : ia=i+np1
555 0 : do 323 nu=ka1,ik1
556 0 : 323 wd=wd+p(ia+(nu-ka-1)*ii)*gd(nu,in1)
557 0 : wa= abs(wd)
558 0 : wb=y(v(j),k)-ucy*wd
559 0 : y(v(j),k)=wb
560 0 : 325 wb= abs(wb)
561 0 : ea(v(j),1)=ea(v(j),1)+wa
562 0 : eb(j)=eb(j)+wb
563 0 : if(wb.lt.1.0e-10) go to 326
564 0 : wb=wa/wb
565 0 : 326 if(wb.lt.ea(v(j),2)) go to 327
566 0 : ea(v(j),2)=wb
567 0 : ea(v(j),3)= float(k)
568 : 327 continue
569 0 : if(i.le.ika) gd(j,in)=wd
570 0 : 328 continue
571 : c
572 : c reset storage indices
573 0 : i=in
574 0 : in=in1
575 0 : in1=i
576 0 : 330 continue
577 : c
578 : c end of second outer loop
579 : c
580 : c iteration complete
581 : c **********
582 : c
583 : c
584 : c set mean relative corrections
585 0 : 340 do 341 i=1,ii
586 0 : 341 ea(v(i),1)=ea(v(i),1)/eb(i)
587 : c set det
588 0 : if(detsgn) 345,350,350
589 0 : 345 det=det*1.e10
590 : c
591 : c
592 0 : 350 return
593 : c
594 : c
595 : c diagnostics
596 : c **********
597 0 : 400 write(istdou,1001)
598 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdou,1001)
599 0 : 401 do 402 k=1,2
600 0 : do 402 j=1,ik
601 0 : write(istdou,1002) (gs(i,j,k), i=1,kap)
602 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdpr,1002)
603 0 : * (gs(i,j,k), i=1,kap)
604 0 : 402 continue
605 0 : write(istdou,1101) (v(i), i=1,ii)
606 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdpr,1101)
607 0 : * (v(i), i=1,ii)
608 0 : go to 500
609 0 : 403 write(istdou,1003)n
610 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdpr,1003)n
611 0 : go to 500
612 0 : 404 write(istdou,1004)
613 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdpr,1004)
614 0 : do 405 i=1,kap
615 0 : do 405 j=1,ik
616 0 : do 405 k=1,2
617 0 : 405 gs(i,j,k)=0.0
618 0 : call bc(x(1),x(nn),y,y(1,nn),zk,ap,aq,g,gs,ik,ik1,nn)
619 0 : go to 401
620 0 : 407 write(istdou,1006) nn,x(1)
621 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdou,1006) nn,x(1)
622 0 : go to 500
623 0 : 408 dx=2.0*dx
624 0 : write(istdou,1007) n,x(n),np1,x(np1),dx,r,nn,(x(i), i=1,nn)
625 0 : if(istdpr.ne.istdou.and.istdpr.gt.0) write(istdpr,1007)
626 0 : * n,x(n),np1,x(np1),dx,r,nn,(x(i), i=1,nn)
627 : c
628 0 : 500 v(1)=0
629 0 : return
630 : c
631 : c
632 : 1000 format(//1x,10('*'),10x,'improper formulation detected by nrk',
633 : . 10x,10('*')/17x,'ii =',i3,4x,'kk =',i3,4x,'ka =',i3,4x,'kb =',i3,
634 : . 4x,'ki =',i3/)
635 : 1001 format(//1x,10('*'),5x,'first boundary condition matrix singular i
636 : .n nrk',5x,10('*')//' gd ='/)
637 : 1002 format(1x,1p14e9.1/)
638 : 1003 format(//1x,10('*'),5x,'e matrix singular at n =',i4,' in nrk',
639 : . 5x,10('*'))
640 : 1004 format(//1x,10('*'),5x,'second boundary matrix singular in nrk',
641 : . 5x,10('*')//' gd ='/)
642 : 1006 format(//1x,10('*'),5x,'null range of independent variable in ',
643 : . 'nrk',5x,10('*')//21x,'x(1) = x(',i4,') =',1pe14.6/)
644 : 1007 format(//1x,10('*'),5x,'independent variable not monotonic in',
645 : . ' nrk',5x,10('*')//16x, 'x(',i4,') =',1pe12.4,
646 : . ', x(',i4,') =',1pe12.4/16x,'difference =',1pe12.4,
647 : . ', range =',1pe12.4//1x,'x(n), n=1,',i4,/
648 : . (1x,1p10e13.5))
649 : c
650 : c warning message
651 : 1100 format(1x,5('*'),3x,'v set in nrk v(',i2,') =',i3)
652 : 1101 format(/9x,'v =',7(1x,5i3))
653 : c
654 : c
655 : end
656 : c.. subroutine leqsd(a,b,n,m,ia,ib,det,isa,isb)
657 : c.. implicit double precision (a-h,o-z)
658 : c.. dimension a(ia,1),b(ib,1)
659 : c..c
660 : c.. call leq(a,b,n,m,ia,ib,det)
661 : c.. return
662 : c.. end
|
