← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/poly/wdmpmu.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/wdmpmu.f
 
1
C/MEMBR ADD NAME=WDMPMU,SSI=0
 
2
c     Copyright INRIA
 
3
      subroutine wdmpmu(mp1r,mp1i,d1,nl1,mp2r,d2,nl2,
 
4
     & mp3r,mp3i,d3,l,m,n)
 
5
c!purpose
 
6
c  ce sous programme effectue le calcul du produit de matrices
 
7
c de polynomes a coefficients complexes
 
8
c
 
9
c           mp3 = mp1 * mp2
 
10
c
 
11
c! calling sequence
 
12
c
 
13
c     pm1 : tableau contenant les coefficients des polynomes,
 
14
c           le coefficient de degre k du polynome pm1(i,j) est range
 
15
c           dans pm1( d1(i + (j-1)*nl1 + k) )
 
16
c           pm1 doit etre de taille au moins d1(nl1*n+1)-d1(1)
 
17
c     d1 : tableau entier de taille nl1*n+1,  si k=i+(j-1)*nl1 alors
 
18
c          d1(k)) contient  l'adresse dans pm1 du coeff de degre 0
 
19
c          du polynome pm1(i,j). Le degre du polynome pm1(i,j) vaut:
 
20
c          d1(k+1)-d1(k) -1
 
21
c     nl1 : entier definissant le rangement dans d1
 
22
c
 
23
c     pm2,d2,nl2 : definitions similaires a celles de pm1,d1,nl1
 
24
c     pm3,d3 : definitions similaires a celles de pm1 et d1, nl3 est
 
25
c              suppose egal a l
 
26
c     l : nombre de lignes de pm1
 
27
c     m : nombre de ligne des matrices pm
 
28
c     n : nombre de colonnes des matrices pm
 
29
c avec les conventions suivantes:
 
30
c
 
31
c     -si l =  0 signifie  que  mp1 est un polynome, et que la
 
32
c multiplication s'entend  alors  au sens de la multiplication
 
33
c de tous les coefficients de mp2 par mp1.
 
34
c
 
35
c     -si n =  0 signifie que mp2 est un polynome, et  que  la
 
36
c multiplication s'entend au sens de la multiplication de tous
 
37
c les coefficients de mp1 par mp2.
 
38
c
 
39
c     -si m =  0 la multiplication s'entend element par element,
 
40
c il est alors suppose que mp1 et mp2, donc mp3 sont de dimension
 
41
c l x n.
 
42
c
 
43
c!
 
44
c auteur: c. klimann, inria
 
45
c date: 25 fevrier 1985.
 
46
c
 
47
c
 
48
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 
49
c
 
50
c
 
51
      double precision mp1r(*),mp1i(*),mp2r(*),mp3r(*),mp3i(*)
 
52
      integer d1(*), d2(*), d3(*)
 
53
      integer l, m, n
 
54
      integer nl1, nl2
 
55
      integer i, j, k, k1, k2, k3
 
56
      integer p1, p2, p3
 
57
c
 
58
c
 
59
      d3(1)= 1
 
60
      if (l .eq. 0 .or. m .eq. 0 .or. n .eq. 0) goto 500
 
61
c
 
62
      p2 = -nl2
 
63
      p3 = -l
 
64
      do 10 j = 1, n
 
65
      p2 = p2 + nl2
 
66
      p3 = p3 + l
 
67
         do 11 i = 1, l
 
68
            mp3r(d3(p3+i)) = 0.0d+0
 
69
            mp3i(d3(p3+i)) = 0.0d+0
 
70
            k3 = 0
 
71
            p1 = i - nl1
 
72
            do 12 k = 1, m
 
73
               p1 = p1 + nl1
 
74
               k2 = d2(p2+k+1) - d2(p2+k) - 1
 
75
               k1 = d1(p1 + 1) - d1(p1) - 1
 
76
               kk=k3
 
77
               call dpmul(mp1r(d1(p1)),k1,mp2r(d2(p2+k)),k2,
 
78
     &         mp3r(d3(p3+i)),kk)
 
79
               call dpmul(mp1i(d1(p1)),k1,mp2r(d2(p2+k)),k2,
 
80
     &         mp3i(d3(p3+i)),k3)
 
81
   12       continue
 
82
            d3(p3 + i + 1) = d3(p3 + i) + k3 + 1
 
83
   11 continue
 
84
   10 continue
 
85
      return
 
86
  500 if (l .eq. 0) goto 600
 
87
      if (m .eq. 0) goto 700
 
88
      p1 = -nl1
 
89
      p3 = -l
 
90
      k2 = d2(2) - d2(1) - 1
 
91
      do 510 j = 1, m
 
92
      p1 = p1 + nl1
 
93
      p3 = p3 + l
 
94
         do 510 i = 1, l
 
95
            k3 = 0
 
96
            k1 = d1(p1 + i + 1) - d1(p1 + i) - 1
 
97
            mp3r(d3(p3 + i)) = 0.0d+0
 
98
            kk=k3
 
99
            call dpmul(mp1r(d1(p1+i)),k1,mp2r(1),k2,mp3r(d3(p3+i)),kk)
 
100
            mp3i(d3(p3 + i)) = 0.0d+0
 
101
            call dpmul(mp1i(d1(p1+i)),k1,mp2r(1),k2,mp3i(d3(p3+i)),k3)
 
102
            d3(p3 + i + 1) = d3(p3 + i) + k3 + 1
 
103
510   continue
 
104
      return
 
105
600   k1 = d1(2) - d1(1) - 1
 
106
      p2 = -nl2
 
107
      p3 = -m
 
108
      do 610 j = 1, n
 
109
      p2 = p2 + nl2
 
110
      p3 = p3 + m
 
111
         do 610  i = 1, m
 
112
            k3 = 0
 
113
            k2 = d2(p2 + i + 1) - d2(p2 + i) - 1
 
114
            mp3r(d3(p3+i)) = 0.0d+0
 
115
            kk=k3
 
116
            call dpmul(mp1r(1),k1,mp2r(d2(p2+i)),k2,mp3r(d3(p3+i)),kk)
 
117
            mp3i(d3(p3+i)) = 0.0d+0
 
118
            call dpmul(mp1i(1),k1,mp2r(d2(p2+i)),k2,mp3i(d3(p3+i)),k3)
 
119
            d3(p3 + i + 1) = d3(p3 + i) + k3 + 1
 
120
610   continue
 
121
      return
 
122
  700 p1 = -nl1
 
123
      p2 = -nl2
 
124
      p3 = -l
 
125
      do 710 j = 1, n
 
126
         p1 = p1 + nl1
 
127
         p2 = p2 + nl2
 
128
         p3 = p3 + l
 
129
         do 710 i = 1, l
 
130
            k1 = d1(p1 + i + 1) - d1(p1 + i) - 1
 
131
            k2 = d2(p2 + i + 1) - d2(p2 + i) - 1
 
132
            mp3r(d3(p3+i)) = 0.0d+0
 
133
            k3 = 0
 
134
            call dpmul(mp1r(d1(p1+i)),k1,mp2r(d2(p2+i)),k2,
 
135
     1                 mp3r(d3(p3+i)),k3)
 
136
            mp3i(d3(p3+i)) = 0.0d+0
 
137
            k3 = 0
 
138
            call dpmul(mp1i(d1(p1+i)),k1,mp2r(d2(p2+i)),k2,
 
139
     1                 mp3i(d3(p3+i)),k3)
 
140
            d3(p3 + i + 1) = d3(p3 + i) + k3  + 1
 
141
710   continue
 
142
      return
 
143
      end