← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/poly/dpmul.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/dpmul.f
 
1
C/MEMBR ADD NAME=DPMUL,SSI=0
 
2
c     Copyright INRIA
 
3
      subroutine dpmul(p1,d1,p2,d2,p3,d3)
 
4
c!purpose
 
5
c  ce sous programme effectue le calcul:
 
6
c
 
7
c                p3(x) = p3(x) + (p1(x) * p2(x))
 
8
c
 
9
c ou  p1,  p2  et p3 sont des polynomes de degres d1, d2 et d3,
 
10
c respectivement. Ils sont classes en ordre croissant.
 
11
c Tous  les  parametres  sont  d'entree, sauf p3 et d3 qui sont
 
12
c d'entree-sortie.
 
13
c!
 
14
c auteur: c. klimann.
 
15
c date: 22 fevrier 1985.
 
16
c &&var
 
17
      implicit double precision (a-h,o-z)
 
18
      double precision dlamch,ddot
 
19
      double precision p1(*),p2(*),p3(*)
 
20
      integer d1,d2,d3
 
21
      integer dmax,dmin,dsum
 
22
      integer i, k, j, l
 
23
      integer e1, e2
 
24
      eps=dlamch('p')
 
25
c &&ker
 
26
      dsum = d1 + d2
 
27
c fixation de dmax et dmin
 
28
      dmax = d1
 
29
      if (d2 .gt. d1) dmax = d2
 
30
      dmin = dsum - dmax
 
31
c fixation de d3
 
32
      if (d3 .ge. dsum) goto 1
 
33
      e1 = d3+2
 
34
      e2 = dsum +1
 
35
      do 2 i = e1 , e2
 
36
         p3(i) = 0.0d+0
 
37
    2 continue
 
38
      d3 = dsum
 
39
    1 continue
 
40
c cas des degres egaux a zero
 
41
      if (d1 .eq. 0 .or. d2 .eq. 0) goto 53
 
42
c produit et addition
 
43
      e1 = 1
 
44
      e2 = dmin + 1
 
45
      do 10 i = e1, e2
 
46
         w=ddot(i,p1(1),1,p2(1),-1)
 
47
         w1=p3(i)+w
 
48
         if(abs(w1).gt.eps*max(abs(p3(i)),abs(w))) then
 
49
            p3(i) = w1
 
50
         else
 
51
            p3(i)=0.0d+0
 
52
         endif
 
53
   10 continue
 
54
      k = 1
 
55
      if (d1 .eq. d2) goto 21
 
56
      e1 = dmin + 2
 
57
      e2 = dmax + 1
 
58
c si d1 > d2
 
59
      if (d1 .lt. d2) goto 25
 
60
      do 20 i = e1, e2
 
61
        k = k + 1
 
62
        w=ddot(dmin+1, p1(k), 1, p2(1), -1)
 
63
        w1=p3(i)+w
 
64
        if(abs(w1).gt.eps*max(abs(p3(i)),abs(w))) then
 
65
           p3(i) = w1
 
66
        else
 
67
           p3(i)=0.0d+0
 
68
        endif
 
69
   20 continue
 
70
   21 e1 = dmax + 2
 
71
      e2 = dsum + 1
 
72
      l = 1
 
73
      j = dmin + 1
 
74
      do 40 i = e1, e2
 
75
        j = j -1
 
76
        k = k + 1
 
77
        l =l + 1
 
78
        w=ddot(j, p1(k), 1, p2(l), -1)
 
79
        w1=p3(i)+w
 
80
        if(abs(w1).gt.eps*max(abs(p3(i)),abs(w))) then
 
81
           p3(i) = w1
 
82
        else
 
83
           p3(i)=0.0d+0
 
84
        endif
 
85
   40 continue
 
86
      return
 
87
c si d2 > d1
 
88
   25 do 30 i = e1, e2
 
89
        k = k + 1
 
90
        w=ddot(dmin+1, p2(k), -1, p1(1), 1)
 
91
        w1=p3(i)+w
 
92
        if(abs(w1).gt.eps*max(abs(p3(i)),abs(w))) then
 
93
           p3(i) = w1
 
94
        else
 
95
           p3(i)=0.0d+0
 
96
        endif
 
97
   30 continue
 
98
      e1 = dmax + 2
 
99
      e2 = dsum + 1
 
100
      l = 1
 
101
      j = dmin + 1
 
102
      do 50 i = e1, e2
 
103
        j = j -1
 
104
        k = k + 1
 
105
        l =l + 1
 
106
        w=ddot(j, p1(l), 1, p2(k), -1)
 
107
        w1=p3(i)+w
 
108
        if(abs(w1).gt.eps*max(abs(p3(i)),abs(w))) then
 
109
           p3(i) = w1
 
110
        else
 
111
           p3(i)=0.0d+0
 
112
        endif
 
113
   50 continue
 
114
      return
 
115
c cas des degres egaux a zero
 
116
   53 if (d1 .eq. 0 .and. d2 .eq. 0) goto 100
 
117
      e1 = 1
 
118
      if (d1 .eq.0) goto 60
 
119
      e2 = d1 + 1
 
120
      do 55 i = e1, e2
 
121
        w=p1(i) * p2(1)
 
122
        w1=p3(i)+w
 
123
        if(abs(w1).gt.eps*max(abs(p3(i)),abs(w))) then
 
124
           p3(i) = w1
 
125
        else
 
126
           p3(i)=0.0d+0
 
127
        endif
 
128
   55 continue
 
129
      return
 
130
   60 e2 = d2 + 1
 
131
      do 65 i = e1, e2
 
132
         w=p2(i) * p1(1)
 
133
        w1=p3(i)+w
 
134
        if(abs(w1).gt.eps*max(abs(p3(i)),abs(w))) then
 
135
           p3(i) = w1
 
136
        else
 
137
           p3(i)=0.0d+0
 
138
        endif
 
139
   65 continue
 
140
      return
 
141
  100 p3(1) = p3(1) + p1(1) * p2(1)
 
142
      return
 
143
      end