← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/control/hessl2.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/control/hessl2.f
 
1
      subroutine hessl2(neq,tq,pd,nrowpd)
 
2
c!but
 
3
c     Elle etablit la valeur de la Hessienne, derivee
 
4
c       seconde de la fonction phi au point q .
 
5
c!liste d'appel
 
6
c     subroutine hessl2(neq,tq,pd,nrowpd)
 
7
c     Entree :
 
8
c     - neq. tableau entier de taille 3+(nq+1)*(nq+2)
 
9
c         neq(1)=nq est le degre effectif du polynome tq (ou q).
 
10
c         neq(2)=ng est le nombre de coefficient de fourier
 
11
c         neq(3)=dgmax degre maximum pour q (l'adresse des coeff de fourier dans
 
12
c               tq est neq(3)+2
 
13
c         neq(4:(nq+1)*(nq+2)) tableau de travail entier
 
14
c     - tq. tableau reel de taille au moins
 
15
c               6+dgmax+5*nq+6*ng+nq*ng+nq**2*(ng+1)
 
16
c         tq(1:nq+1) est le tableau des coefficients du polynome.
 
17
c         tq(dgmax+2:dgmax+ng+2) est le tableau des coefficients
 
18
c                      de fourier
 
19
c         tq(dgmax+ng+3:) est un tableau de travail de taille au moins
 
20
c                         5+5*nq+5*ng+nq*ng+nq**2*(ng+1)
 
21
c     Sortie :
 
22
c     - pd matrice hessienne
 
23
c!
 
24
c     Copyright INRIA
 
25
      implicit double precision (a-h,o-y)
 
26
      dimension tq(*),pd(nrowpd,*)
 
27
      dimension neq(*)
 
28
c
 
29
      nq=neq(1)
 
30
      ng=neq(2)
 
31
c
 
32
c     decoupage du tableau neq
 
33
      jmxnv=4
 
34
      jmxnw=jmxnv+(nq+1)
 
35
      jw=jmxnw+(nq+1)**2
 
36
c
 
37
c     decoupage du tableau tq
 
38
      itq=1
 
39
      itg=itq+neq(3)+1
 
40
      itr=itg+ng+1
 
41
      itp=itr+nq+ng+1
 
42
      itv=itp+nq+ng+1
 
43
      itw=itv+nq+ng+1
 
44
      itij=itw+nq+ng+1
 
45
      id1aux=itij+ng+1
 
46
      id2aux=id1aux+(ng+1)*nq
 
47
      iw=id2aux+nq*nq*(ng+1)
 
48
 
 
49
      call hl2(nq,tq,tq(itg),ng,pd,nrowpd,tq(itr),
 
50
     $     tq(itp),tq(itv),tq(itw),tq(itij),tq(id1aux),tq(id2aux),
 
51
     $     neq(jmxnv),neq(jmxnw))
 
52
      return
 
53
      end
 
54
 
 
55
 
 
56
 
 
57
      subroutine hl2(nq,tq,tg,ng,pd,nrowpd,tr,tp,tv,tw,tij,d1aux,d2aux,
 
58
     &     maxnv,maxnw)
 
59
c!
 
60
 
 
61
      implicit double precision (a-h,o-y)
 
62
      dimension tq(nq+1),tg(ng+1),pd(nrowpd,*)
 
63
c
 
64
      dimension tr(nq+ng+1),tv(nq+ng+1),tp(nq+ng+1),tw(nq+ng+1),
 
65
     &     tij(ng+1),d1aux(ng+1,nq),d2aux(nq,nq,ng+1)
 
66
      integer maxnv(nq),maxnw(nq,nq)
 
67
c
 
68
c     --- Calcul des derivees premieres de 'vq' ---
 
69
c
 
70
      do 20 i=1,nq
 
71
         if (i.eq.1) then
 
72
c     .     division euclidienne de z^nq*g par q
 
73
            call dset(nq,0.0d0,tp,1)
 
74
            call dcopy(ng+1,tg,1,tp(nq+1),1)
 
75
            call dpodiv(tp,tq,nq+ng,nq)
 
76
            nv1=ng
 
77
c     .     calcul de Lq et Vq
 
78
            call lq(nq,tq,tr,tg,ng)
 
79
            ltvq=nq+1
 
80
c     .     division euclidienne de Vq par q
 
81
            call dcopy(ng+1,tr(ltvq),1,tv,1)
 
82
            call dset(nq,0.0d0,tv(ng+2),1)
 
83
            call dpodiv(tv,tq,ng,nq)
 
84
            nv2=ng-nq
 
85
         else
 
86
            ichoi1=1
 
87
            call dzdivq(ichoi1,nv1,tp,nq,tq)
 
88
            ichoi2=2
 
89
            call mzdivq(ichoi2,nv2,tv,nq,tq)
 
90
         endif
 
91
         maxnv(i)=max(nv1,nv2)
 
92
         do 10 j=1,maxnv(i)+1
 
93
            d1aux(j,i)= tp(nq+j)-tv(nq+j)
 
94
 10      continue
 
95
 20   continue
 
96
c
 
97
c     --- Calcul des derivees secondes de 'vq' ---
 
98
c
 
99
      do 50 i=1,nq
 
100
         call dset(ng+nq+1,0.0d0,tw,1)
 
101
         do 40 j=nq,1,-1
 
102
            if (j.eq.nq) then
 
103
               call dcopy(maxnv(i)+1,d1aux(1,i),1,tw(nq),1)
 
104
               nw=maxnv(i)+nq-1
 
105
               call dpodiv(tw,tq,nw,nq)
 
106
               nw=nw-nq
 
107
            else
 
108
               ichoix=1
 
109
               call dzdivq(ichoix,nw,tw,nq,tq)
 
110
            endif
 
111
            do 30 k=1,nw+1
 
112
               d2aux(i,j,k)=tw(nq+k)
 
113
 30         continue
 
114
            maxnw(i,j)=nw
 
115
 40      continue
 
116
 50   continue
 
117
c
 
118
c     --- Conclusion des calculs sur la hessienne ---
 
119
c
 
120
      do 100 i=1,nq
 
121
         do 90 j=1,i
 
122
            call scapol(maxnv(i),d1aux(1,i),maxnv(j),
 
123
     &           d1aux(1,j),y1)
 
124
c     
 
125
            if (maxnw(i,j).gt.maxnw(j,i)) then
 
126
               maxij=maxnw(i,j)
 
127
               minij=maxnw(j,i)
 
128
               do 60 k=minij+2,maxij+1
 
129
                  tij(k)= -d2aux(i,j,k)
 
130
 60            continue
 
131
            else if (maxnw(i,j).lt.maxnw(j,i)) then
 
132
               maxij=maxnw(j,i)
 
133
               minij=maxnw(i,j)
 
134
               do 70 k=minij+2,maxij+1
 
135
                  tij(k)= -d2aux(j,i,k)
 
136
 70            continue
 
137
            else
 
138
               maxij=maxnw(i,j)
 
139
               minij=maxij
 
140
            endif
 
141
c     
 
142
            do 80 k=1,minij+1
 
143
               tij(k)= -d2aux(i,j,k) -d2aux(j,i,k)
 
144
 80         continue
 
145
c     
 
146
            call scapol(maxij,tij,ng,tr(ltvq),y2)
 
147
 
 
148
            if (i.eq.j) then
 
149
               pd(i,i)=-2.0d+0 * (y1+y2)
 
150
            else
 
151
               pd(i,j)=-2.0d+0 * (y1+y2)
 
152
               pd(j,i)=-2.0d+0 * (y1+y2)
 
153
            endif
 
154
 90      continue
 
155
 100  continue
 
156
      return
 
157
      end