← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/poly/dmrdsp.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/poly/dmrdsp.f
 
1
      subroutine dmrdsp(mpn,dn,mpd,dd,nl,mm,nn,var,lvar,maxc,mode,ll,
 
2
     1    lunit,cw,iw)
 
3
c!but 
 
4
c     dmpdsp ecrit une matrice polynomiale (ou un polynome) sous
 
5
c     la forme d'un tableau de polynomes, avec gestion automatique de
 
6
c     l'espace disponible.
 
7
c!liste d'appel
 
8
c     
 
9
c     subroutine dmrdsp(mpn,dn,mpd,dd,nl,m,n,var,lvar,maxc,mode,ll,
 
10
c     1 lunit,cw,iw)
 
11
c     
 
12
c     double precision mp(*)
 
13
c     integer d(nl*n+1),nl,m,n,lvar,maxc,mode,iw(*),ll,lunit
 
14
c     character var*(*),cw*(*)
 
15
c     
 
16
c     pm : tableau reel contenant les coefficients des polynomes,
 
17
c     le coefficient de degre k du polynome pm(i,j) est range
 
18
c     dans pm( d(i + (j-1)*nl + k) )
 
19
c     pm doit etre de taille au moins d(nl*n+1)-d(1)
 
20
c     d : tableau entier de taille nl*n+1,  si k=i+(j-1)*nl alors
 
21
c     d(k)) contient  l'adresse dans pm du coeff de degre 0
 
22
c     du polynome pm(i,j). Le degre du polynome pm(i,j) vaut:
 
23
c     d(k+1)-d(k) -1
 
24
c     nl : entier definissant le rangement dans d
 
25
c     m : nombre de ligne de la matrice polynomiale
 
26
c     n : nombre de colonnes de la matrice polynomiale
 
27
c     var : nom de la variable muette
 
28
c     lvar : nombre de caracteres de var
 
29
c     maxc : nombre de caracteres maximum autorise pour
 
30
c     representer un nombre
 
31
c     mode : si mode = 1 alors representation variable
 
32
c     si mode = 0 alors representation d(maxc).(maxc-7)
 
33
c     ll : longueur de ligne maximum admissible
 
34
c     lunit : etiquette logique du support d'edition
 
35
c     cw : chaine de caracteres de travail de longueur au moins ll*2
 
36
c     iw : tableau de travail entier de taille au moins egale a
 
37
c     n*(4+m)+1+dn(n*m+1)+dd(n*m+1)
 
38
c!origine
 
39
c     s. steer inria 1986
 
40
c     Copyright INRIA
 
41
 
 
42
c!    
 
43
c     
 
44
      double precision mpn(*),mpd(*),a
 
45
      integer dd(*),dn(*),iw(*),maxc,mode
 
46
      integer fl,sk,sl,c1,c2,typ
 
47
      logical first
 
48
      character var*(*),cw*(*),sgn*1,dl*1
 
49
      character*10 form(2),fexp,expo
 
50
      integer nind
 
51
c     
 
52
      data nind/5/
 
53
c     
 
54
      m=abs(mm)
 
55
      n=abs(nn)
 
56
c
 
57
      cw=' '
 
58
      write(form(1),130) maxc,maxc-7
 
59
      dl=' '
 
60
      if(m*n.gt.1) dl='!'
 
61
c     
 
62
c     phase d'analyse: pour chaque coefficient a representer on determine
 
63
c     format avec lequel on va l'editer, on en deduit la longueur
 
64
c     de la representation de chacun des polynomes.
 
65
c     les differents formats sont stockes sous forme codee dans iw
 
66
c     a partir de lf
 
67
c     la taille respective des representation des chacun des polynomes
 
68
c     est contenue dans iw a partir de 1 .
 
69
c     
 
70
      lcol=1
 
71
      lbloc=lcol+n-1
 
72
      lfn=lbloc+n+2
 
73
      lfd=lfn+dn(n*m+1)
 
74
      ldelta=lfd+dd(n*m+1)
 
75
      ldeb=ldelta+m*n
 
76
      lfin=ldeb+n
 
77
c     
 
78
      lines=0
 
79
      nbloc=1
 
80
      iw(lbloc+nbloc)=n
 
81
      sk=0
 
82
      ldefn=lfn
 
83
      ldg=-nl
 
84
      ldefd=lfd
 
85
      idelta=ldelta
 
86
c     
 
87
      k0=1
 
88
      do 11 k=1,n
 
89
         sl=0
 
90
         iw(lcol-1+k)=0
 
91
         ldg=ldg+nl
 
92
         do 10 l=1,m
 
93
c     
 
94
c     traitement du polynome (l,k)
 
95
            lpn=dn(ldg+l)-1
 
96
            npn=dn(ldg+l+1)-dn(ldg+l)
 
97
            lghn=0
 
98
            first=.true.
 
99
            do 05 i=1,npn
 
100
               a=abs(mpn(lpn+i))
 
101
               iw(ldefn)=0
 
102
               if(a.ne.0.0d+0) then
 
103
                  first=.false.
 
104
c     determination du format devant representer a
 
105
                  typ=1
 
106
                  if(mode.eq.1) call fmt(a,maxc,typ,n1,n2)
 
107
                  if(typ.eq.2) then
 
108
                     fl=n1
 
109
                     iw(ldefn)=n2+32*n1
 
110
                  elseif(typ.lt.0) then
 
111
                     iw(ldefn)=typ
 
112
                     n2=1
 
113
                     fl=3
 
114
                  else
 
115
                     iw(ldefn)=1
 
116
                     fl=maxc
 
117
                     n2=maxc-7
 
118
                  endif
 
119
c     
 
120
c     determination de la longueur de la representation du monome,
 
121
c     cette longueur est a priori fl+2 (' '//sgn//rep(a)//var).
 
122
c     mais peut etre reduite dans des cas particulier
 
123
                  lghn=lghn+fl+2
 
124
                  if(n2.eq.0) then
 
125
                     lghn=lghn-1
 
126
                     if(i.ne.1.and.int(a+0.1).eq.1) lghn=lghn-1
 
127
                  endif
 
128
                  if(i.ne.1) lghn=lghn+lvar
 
129
               endif
 
130
               ldefn=ldefn+1
 
131
 05         continue
 
132
            
 
133
c     
 
134
c     cas particulier du dernier exposant du polynome
 
135
            nd=ifix(log10(0.5+npn))+1
 
136
            lghn=lghn+nd
 
137
c     cas particulier d'un polynome reduit a 0
 
138
            if(first) lghn=4
 
139
c     
 
140
            lpd=dd(ldg+l)-1
 
141
            npd=dd(ldg+l+1)-dd(ldg+l)
 
142
            lghd=0
 
143
            first=.true.
 
144
            do 08 i=1,npd
 
145
               a=abs(mpd(lpd+i))
 
146
               iw(ldefd)=0
 
147
               if(a.ne.0.0d+0) then
 
148
                  first=.false.
 
149
c     determination du format devant representer a
 
150
                  typ=1
 
151
                  if(mode.eq.1) call fmt(a,maxc,typ,n1,n2)
 
152
                  if(typ.eq.2) then
 
153
                     fl=n1
 
154
                     iw(ldefd)=n2+32*n1
 
155
                  elseif(typ.lt.0) then
 
156
                     iw(ldefd)=typ
 
157
                     n2=1
 
158
                     fl=3
 
159
                  else
 
160
                     iw(ldefd)=1
 
161
                     fl=maxc
 
162
                     n2=maxc-7
 
163
                  endif
 
164
c     determination de la longueur de la representation du monome,
 
165
c     cette longueur est a priori fl+2 (' '//sgn//rep(a)//var).
 
166
c     mais peut etre reduite dans des cas particulier
 
167
                  lghd=lghd+fl+2
 
168
                  if(n2.eq.0) then
 
169
                     lghd=lghd-1
 
170
                     if(i.ne.1.and.int(a+0.1).eq.1) lghd=lghd-1
 
171
                  endif
 
172
                  if(i.ne.1) lghd=lghd+lvar
 
173
               endif
 
174
               ldefd=ldefd+1
 
175
 08         continue
 
176
c     
 
177
c     cas particulier du dernier exposant du polynome
 
178
            nd=ifix(log10(0.5+npd))+1
 
179
            lghd=lghd+nd
 
180
c     cas particulier d'un polynome reduit a 0
 
181
            if(first) lghd=4
 
182
c     
 
183
            iw(k)=max(iw(k),lghn,lghd)
 
184
            sl=sl+(lghn/(ll-10))+(lghd/(ll-10))+2
 
185
            iw(idelta)=min(lghn,ll-2)-min(lghd,ll-2)
 
186
            idelta=idelta+1
 
187
c     
 
188
 10      continue
 
189
         sk=sk+iw(k)
 
190
         if(sk.gt.ll-2) then
 
191
            if(k.eq.k0) then
 
192
               iw(lbloc+nbloc)=k
 
193
               sk=0
 
194
               k0=k+1
 
195
            else
 
196
               iw(lbloc+nbloc)=k-1
 
197
               sk=iw(k)
 
198
               k0=k
 
199
            endif
 
200
            nbloc=nbloc+1
 
201
            iw(lbloc+nbloc)=n
 
202
c     lines=lines+2*sl+m+2
 
203
         endif
 
204
 11   continue
 
205
      nbloc=min(nbloc,n)
 
206
c     
 
207
      l1=1
 
208
      if(mm.lt.0) then
 
209
         write(cw(l1:l1+4),'(''eye *'')') 
 
210
         l1=l1+5
 
211
         call basout(io,lunit,cw(1:l1-1))
 
212
         call basout(io,lunit,' ')
 
213
         if(io.eq.-1) goto 99
 
214
      endif
 
215
 
 
216
c     
 
217
c     phase d'edition : les deux chaines de caracteres representant
 
218
c     la ligne des exposants et la ligne des coefficients,sont
 
219
c     constituees puis imprimees.
 
220
c     
 
221
      k1=1
 
222
      do 70 ib=1,nbloc
 
223
         k2=iw(lbloc+ib)
 
224
         ll1=0
 
225
         if(nbloc.ne.1) then
 
226
            call blktit(lunit,k1,k2,io)
 
227
            if (io.eq.-1) goto 99
 
228
         endif
 
229
c     
 
230
         cw(1:1)=dl
 
231
         c1=2
 
232
         cw(1+ll:1+ll)=dl
 
233
         c2=max(3+ll,nind+maxc+15)
 
234
c     
 
235
         do 60 l=1,m
 
236
c     numerateur
 
237
            l1=c1
 
238
            l2=c2
 
239
            if(iw(k1).gt.ll-2) ll1=ll
 
240
            do 45 k=k1,k2
 
241
               l0=l1
 
242
               idelta=ldelta-1+l+(k-1)*m
 
243
               ndelta=0
 
244
               if(iw(idelta).lt.-1) then
 
245
                  ndelta=-iw(idelta)/2
 
246
                  cw(l1:l1+ndelta-1)=' '
 
247
                  cw(l2:l2+ndelta-1)=' '
 
248
                  l1=l1+ndelta
 
249
                  l2=l2+ndelta
 
250
               endif
 
251
c     
 
252
               ldg=(k-1)*nl+l
 
253
               lpn=dn(ldg)-1
 
254
               npn=dn(ldg+1)-dn(ldg)
 
255
               ldefn=lfn-1+dn(ldg)-dn(1)
 
256
               first=.true.
 
257
c     
 
258
               iw(ldeb-1+k)=l2
 
259
               iw(lfin-1+k)=0
 
260
               do 40 j=1,npn
 
261
                  ifmt=iw(ldefn+j)
 
262
                  if(ifmt.eq.0) goto 40
 
263
                  sgn='+'
 
264
                  if(first) sgn=' '
 
265
                  first=.false.
 
266
                  if(mpn(lpn+j).lt.0.0d+0) sgn='-'
 
267
                  a=abs(mpn(lpn+j))
 
268
c     
 
269
                  if(ifmt.eq.1) then
 
270
                     nf=1
 
271
                     fl=maxc
 
272
                     n2=1
 
273
                  elseif(ifmt.ge.0) then
 
274
                     nf=2
 
275
                     n1=ifmt/32
 
276
                     n2=ifmt-32*n1
 
277
                     fl=n1
 
278
                     write(form(nf),120) fl,n2
 
279
                  elseif(ifmt.lt.0) then
 
280
c     Inf/Nan
 
281
                     fl=3
 
282
                     n2=1
 
283
                  endif
 
284
 
 
285
c     
 
286
                  nd=0
 
287
                  if(j.gt.2) nd=ifix(log10(0.5+j))+1
 
288
                  if(l2+fl+2+lvar+nd.gt.c2+ll-2) then
 
289
c     gestion des lignes suites
 
290
                     if(l1.le.ll-1) then
 
291
                        cw(l1:ll-1)=' '
 
292
                        l1=ll
 
293
                     endif
 
294
                     if(l2.le.c2+ll-2) then
 
295
                        cw(l2:c2+ll-2)=' '
 
296
                        l2=c2+ll-2
 
297
                     endif
 
298
                     iw(lfin-1+k)=l2-1
 
299
                     cw(l1:l1)=dl
 
300
                     call basout(io,lunit,cw(c1-1:l1))
 
301
                     cw(l2:l2)=dl
 
302
                     cw(c2-1:c2-1)=dl
 
303
                     call basout(io,lunit,cw(c2-1:l2))
 
304
                     if(io.eq.-1) goto 99
 
305
                     cw(c2-1:c2+nind-1)=' '
 
306
                     cw(c2-1:c2-1)=dl
 
307
                     l2=c2+nind
 
308
                     cw(c1-1:c1+nind-1)=' '
 
309
                     cw(c1-1:c1-1)=dl
 
310
                     l1=c1+nind
 
311
                  endif
 
312
c     representation du monome
 
313
                  cw(l2:l2+1)=' '//sgn
 
314
                  l2=l2+1
 
315
                  if(ifmt.ge.0) then
 
316
                     write(cw(l2+1:l2+fl),form(nf)) a
 
317
                  elseif(ifmt.eq.-1) then
 
318
                     cw(l2+1:l2+fl)='Inf'
 
319
                  elseif(ifmt.eq.-2) then
 
320
                     cw(l2+1:l2+fl)='Nan'
 
321
                  endif
 
322
                  l2=l2+fl
 
323
                  if(n2.eq.0) l2=l2-1
 
324
                  if(j.gt.1) then
 
325
                     if(n2.eq.0.and.int(a+0.1).eq.1) l2=l2-1
 
326
                     cw(l2+1:l2+lvar)=var(1:lvar)
 
327
                     l2=l2+lvar
 
328
                  endif
 
329
                  nl1=l2+c1-c2
 
330
                  cw(l1:nl1)=' '
 
331
                  if(j.gt.2) then
 
332
                     write(fexp,110) nd
 
333
                     write(expo,fexp) j-1
 
334
                     cw(nl1+1:nl1+nd)=expo(1:nd)
 
335
                     l1=nl1+nd
 
336
                  endif
 
337
                  l1=l1+1
 
338
                  l2=l2+1
 
339
 40            continue
 
340
               if(first) then
 
341
c     cas particulier du polynome nul
 
342
                  cw(l1:l1+3)=' '
 
343
                  cw(l2:l2+3)='   0'
 
344
                  l1=l1+4
 
345
                  l2=l2+4
 
346
                  nd=0
 
347
               endif
 
348
               if(iw(lfin-1+k).eq.0) iw(lfin-1+k)=l2
 
349
               if(nd.ne.0) cw(l2:l2+nd-1)=' '
 
350
               nl1=l0+iw(k)
 
351
               if(ll1.eq.ll) nl1=ll-1
 
352
               cw(l1:nl1)=' '
 
353
               l1=nl1+1
 
354
               cw(l2:c2+nl1-c1)=' '
 
355
               l2=c2+nl1-c1+1
 
356
 45         continue
 
357
            if(cw(c1:l1-1).ne.' ') then
 
358
               cw(l1:l1)=dl
 
359
               call basout(io,lunit,cw(c1-1:l1))
 
360
            endif
 
361
            cw(l2:l2)=dl
 
362
            cw(c2-1:c2-1)=dl
 
363
            call basout(io,lunit,cw(c2-1:l2))
 
364
            if(io.eq.-1) goto 99
 
365
c     
 
366
c     trait de fraction
 
367
            cw(c2:l2-1)=' '
 
368
            jjb1=c2
 
369
            do 47 k=k1,k2
 
370
               idelta=ldelta-1+l+(k-1)*m
 
371
               ndelta=max(0,-iw(idelta)/2)
 
372
               ideb=max(jjb1,iw(ldeb-1+k)-ndelta+2)
 
373
               ifin=iw(lfin-1+k)+ndelta-2
 
374
               if(ifin-ideb+1.eq.2) ideb=ideb-1
 
375
               do 46 i=ideb,ifin
 
376
                  cw(i+1:i+1)='-'
 
377
 46            continue
 
378
               jjb1=iw(lfin-1+k)+1
 
379
 47         continue
 
380
            cw(l2:l2)=dl
 
381
            call basout(io,lunit,cw(c2-1:l2))
 
382
            if(io.eq.-1) goto 99
 
383
c     
 
384
c     denominateur
 
385
            l1=c1
 
386
            l2=c2
 
387
            do 55 k=k1,k2
 
388
               l0=l1
 
389
               idelta=ldelta-1+l+(k-1)*m
 
390
               ndelta=0
 
391
               if(iw(idelta).gt.1) then
 
392
                  ndelta=iw(idelta)/2
 
393
                  cw(l1:l1+ndelta-1)=' '
 
394
                  cw(l2:l2+ndelta-1)=' '
 
395
                  l1=l1+ndelta
 
396
                  l2=l2+ndelta
 
397
               endif
 
398
c     
 
399
               ldg=(k-1)*nl+l
 
400
               lpd=dd(ldg)-1
 
401
               npd=dd(ldg+1)-dd(ldg)
 
402
               ldefd=lfd-1+dd(ldg)-dd(1)
 
403
               first=.true.
 
404
c     
 
405
               do 50 j=1,npd
 
406
                  ifmt=iw(ldefd+j)
 
407
                  if(ifmt.eq.0) goto 50
 
408
                  sgn='+'
 
409
                  if(first) sgn=' '
 
410
                  first=.false.
 
411
                  if(mpd(lpd+j).lt.0.0d+0) sgn='-'
 
412
                  a=abs(mpd(lpd+j))
 
413
c     
 
414
                  if(ifmt.eq.1) then
 
415
                     nf=1
 
416
                     fl=maxc
 
417
                     n2=1
 
418
                  elseif(ifmt.ge.0) then
 
419
                     nf=2
 
420
                     n1=ifmt/32
 
421
                     n2=ifmt-32*n1
 
422
                     fl=n1
 
423
                     write(form(nf),120) fl,n2
 
424
                  elseif(ifmt.lt.0) then
 
425
c     Inf/Nan
 
426
                     fl=3
 
427
                     n2=1
 
428
                  endif
 
429
c     
 
430
                  nd=0
 
431
                  if(j.gt.2) nd=ifix(log10(0.5+j))+1
 
432
                  if(l2+fl+2+lvar+nd.gt.c2+ll-2) then
 
433
c     gestion des lignes suites
 
434
                     if(l1.le.ll-1) then
 
435
                        cw(l1:ll-1)=' '
 
436
                        l1=ll
 
437
                     endif
 
438
                     if(l2.le.c2+ll-2) then
 
439
                        cw(l2:c2+ll-2)=' '
 
440
                        l2=c2+ll-2
 
441
                     endif
 
442
                     cw(l1:l1)=dl
 
443
                     call basout(io,lunit,cw(c1-1:l1))
 
444
                     cw(l2:l2)=dl
 
445
                     call basout(io,lunit,cw(c2-1:l2))
 
446
                     if(io.eq.-1) goto 99
 
447
                     cw(c2:c2-1+nind)=' '
 
448
                     cw(c2-1:c2-1)=dl
 
449
                     l2=c2+nind
 
450
                     cw(c1:c1-1+nind)=' '
 
451
                     cw(c1-1:c1-1)=dl
 
452
                     l1=c1+nind
 
453
                  endif
 
454
c     representation du monome
 
455
                  cw(l2:l2+1)=' '//sgn
 
456
                  l2=l2+1
 
457
                  if(ifmt.ge.0) then
 
458
                     write(cw(l2+1:l2+fl),form(nf)) a
 
459
                  elseif(ifmt.eq.-1) then
 
460
                     cw(l2+1:l2+fl)='Inf'
 
461
                  elseif(ifmt.eq.-2) then
 
462
                     cw(l2+1:l2+fl)='Nan'
 
463
                  endif
 
464
                  l2=l2+fl
 
465
                  if(n2.eq.0) l2=l2-1
 
466
                  if(j.gt.1) then
 
467
                     if(n2.eq.0.and.int(a+0.1).eq.1) l2=l2-1
 
468
                     cw(l2+1:l2+lvar)=var(1:lvar)
 
469
                     l2=l2+lvar
 
470
                  endif
 
471
                  nl1=l2+c1-c2
 
472
                  cw(l1:nl1)=' '
 
473
                  if(j.gt.2) then
 
474
                     write(fexp,110) nd
 
475
                     write(expo,fexp) j-1
 
476
                     cw(nl1+1:nl1+nd)=expo(1:nd)
 
477
                     l1=nl1+nd
 
478
                  endif
 
479
                  l1=l1+1
 
480
                  l2=l2+1
 
481
 50            continue
 
482
               if(first) then
 
483
c     cas particulier du polynome nul
 
484
                  cw(l1:l1+3)=' '
 
485
                  cw(l2:l2+3)='   0'
 
486
                  l1=l1+4
 
487
                  l2=l2+4
 
488
                  nd=0
 
489
               endif
 
490
               if(nd.ne.0) cw(l2:l2+nd-1)=' '
 
491
               nl1=l0+iw(k)
 
492
               if(ll1.eq.ll) nl1=ll-1
 
493
               cw(l1:nl1)=' '
 
494
               l1=nl1+1
 
495
               cw(l2:c2+nl1-c1)=' '
 
496
               l2=c2+nl1-c1+1
 
497
 55         continue
 
498
            if(cw(c1:l1-1).ne.' ') then
 
499
               cw(l1:l1)=dl
 
500
               call basout(io,lunit,cw(c1-1:l1))
 
501
            endif
 
502
            cw(l2:l2)=dl
 
503
            cw(c2-1:c2-1)=dl
 
504
            call basout(io,lunit,cw(c2-1:l2))
 
505
            cw(c1:l1-1)=' '
 
506
            cw(l1:l1)=dl
 
507
            if(io.eq.-1) goto 99
 
508
            if(l.ne.m) then
 
509
               call basout(io,lunit,cw(c1-1:l1))
 
510
               if(io.eq.-1) goto 99
 
511
            endif
 
512
 60      continue
 
513
c     
 
514
         k1=k2+1
 
515
 70   continue
 
516
c     
 
517
 99   return
 
518
c     
 
519
 110  format('(i',i2,')')
 
520
 120  format('(f',i2,'.',i2,')')
 
521
 130  format('(1pd',i2,'.',i2,')')
 
522
c     
 
523
      end