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