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
|
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
c Copyright (C) ????-2008 - INRIA - Serge STEER
c
c This file must be used under the terms of the CeCILL.
c This source file is licensed as described in the file COPYING, which
c you should have received as part of this distribution. The terms
c are also available at
c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
C/MEMBR ADD NAME=WMPAD,SSI=0
c Copyright INRIA
subroutine wmpad(pm1r,pm1i,d1,l1,pm2r,pm2i,d2,l2,pm3r,pm3i,d3,
& m,n)
c!but
c cette subroutine ajoute deux matrices dont les coefficients
c sont des polynomes a coefficients complexes : pm3=pm1+pm2
c!liste d'appel
c
c subroutine wmpad(pm1r,pm1i,d1,l1,pm2r,pm2i,d2,l2,pm3r,pm3i,d3,
c & m,n)
c double precision pm1r(*),pm1i(*),pm2r(*),pm2i(*),pm3r(*),pm3i(*)
c integer d1(l1*n+1),d2(l2*n+1),d3(m*n+1),m,n,l1,l2
c
c pm1 : tableau reel contenant les coefficients des polynomes,
c le coefficient de degre k du polynome pm1(i,j) est range
c dans pm1( d1(i + (j-1)*l1 + k) )
c pm1 doit etre de taille au moins d1(l1*n+1)-d1(1)
c d1 : tableau entier de taille l1*n+1, si k=i+(j-1)*l1 alors
c d1(k)) contient l'adresse dans pm1 du coeff de degre 0
c du polynome pm1(i,j). Le degre du polynome pm1(i,j) vaut:
c d1(k+1)-d1(k) -1
c l1 : entier definissant le rangement dans d1
c
c pm2,d2,l2 : definitions similaires a celles de pm1,d1,l1
c pm3,d3 : definitions similaires a celles de pm1 et d1, l3 est
c suppose egal a m
c m : nombre de ligne des matrices pm
c n : nombre de colonnes des matrices pm
c!
double precision pm1r(*),pm1i(*),pm2r(*),pm2i(*),pm3r(*),pm3i(*)
integer d1(*),d2(*),d3(*),m,n,l1,l2
c
integer n1,n2,n3,mn,i,k
c
mn=m*n
c
d3(1)=1
i1=-l1
i2=-l2
k3=0
c boucle sur les polynomes
do 20 j=1,n
i1=i1+l1
i2=i2+l2
do 20 i=1,m
k1=d1(i1+i)-1
k2=d2(i2+i)-1
n1=d1(i1+i+1)-d1(i1+i)
n2=d2(i2+i+1)-d2(i2+i)
if(n1.gt.n2) goto 15
c
c n1.le.n2
c
do 12 k=1,n1
pm3r(k3+k)=pm1r(k1+k)+pm2r(k2+k)
pm3i(k3+k)=pm1i(k1+k)+pm2i(k2+k)
12 continue
if(n1.eq.n2) goto 14
n3=n1+1
do 13 k=n3,n2
pm3r(k3+k)=+pm2r(k2+k)
13 pm3i(k3+k)=pm2i(k2+k)
14 n3=n2
d3(i+1+(j-1)*m)=d3(i+(j-1)*m)+n3
goto 18
c
c n1.gt.n2
c
15 do 16 k=1,n2
pm3r(k3+k)=pm1r(k1+k)+pm2r(k2+k)
16 pm3i(k3+k)=pm1i(k1+k)+pm2i(k2+k)
n3=n2+1
do 17 k=n3,n1
pm3r(k3+k)=pm1r(k1+k)
17 pm3i(k3+k)=pm1i(k1+k)
n3=n1
d3(i+1+(j-1)*m)=d3(i+(j-1)*m)+n3
c
18 k1=k1+n1
k2=k2+n2
k3=k3+n3
20 continue
return
end
|
