1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine subspace(d,aa,bb,cc,alpham,betam,nev,xini,
20     &  cd,cv,time,rwork,lrw,m,jout,rpar,bj,iwork,liw,iddebdf,bjp)
21!
22!     solves the linear dynamic equations mapped on the subspace
23!     of the eigenvectors (only if there are dashpots in the
24!     model)
25!
26      implicit none
27!
28      save time0
29!
30      integer nev,nev2,info(15),idid,lrw,iwork(*),liw,jout,id,
31     &  iab,iaa,ibb,i,j,m,iddebdf
32!
33      real*8 d(*),aa(nev),bb(nev),cc(nev,*),alpham,betam,
34     &  xini(*),cd(*),cv(*),time,time0,rtol,atol,rwork(*),rpar(*),
35     &  bj(*),bjp(*)
36!
37      external df,djac
38!
39!
40      nev2=2*nev
41!
42!     transferring fields into global field rpar
43!     (needed for subroutine fd)
44!     rpar contains (field, size): m+0.5, 1
45!                                  alpham, 1
46!                                  betam, 1
47!                                  cc, nev**2
48!                                  d, nev
49!                                  time
50!                                  aa(1)..aa(nev), nev
51!                                  bb(1)..bb(nev), nev
52!
53      if(m.eq.1) then
54         rpar(2)=alpham
55         rpar(3)=betam
56         do j=1,nev
57            do i=1,nev
58               rpar(3+(j-1)*nev+i)=cc(i,j)
59            enddo
60         enddo
61         id=3+nev*nev
62         do i=1,nev
63            rpar(id+i)=d(i)
64         enddo
65!
66!        copying the initial conditions for the system of first order
67!        differential equations
68!
69         do i=1,nev
70            xini(i)=cd(i)
71            xini(nev+i)=cv(i)
72         enddo
73       endif
74!
75      iaa=3+nev*(1+nev)+1
76      rpar(iaa)=time
77      ibb=iaa+nev
78      do i=1,nev
79         rpar(iaa+i)=aa(i)
80         rpar(ibb+i)=bb(i)
81      enddo
82!
83      do i=1,3
84         info(i)=0
85      enddo
86      info(4)=1
87      info(5)=1
88      info(6)=0
89      rwork(1)=time
90!
91!     absolute and relative tolerance for dderkf
92!
93      rtol=1.d-5
94      atol=1.d-3
95!
96      if(iddebdf.eq.0) then
97         call ddeabm(df,nev2,time0,xini,time,info,rtol,atol,idid,rwork,
98     &        lrw,iwork,liw,rpar,nev)
99!
100         if((idid.ne.2).and.(idid.ne.3)) then
101            write(*,*)
102     &         '*WARNING in subspace: ddeabm did not converge properly'
103            write(*,*) '         idid= ',idid
104            write(*,*) '         switch to routine ddebdf'
105            iddebdf=2
106            return
107         endif
108      else
109         call ddebdf(df,nev2,time0,xini,time,info,rtol,atol,idid,rwork,
110     &        lrw,iwork,liw,rpar,nev,djac)
111         if((idid.ne.2).and.(idid.ne.3)) then
112            write(*,*)
113     &           '*ERROR in subspace: ddebdf did not converge properly'
114            write(*,*) '       idid= ',idid
115            call exit(201)
116         endif
117      endif
118!
119!     copying the solution into field bj
120!
121      do i=1,nev
122         bj(i)=xini(i)
123         bjp(i)=xini(nev+i)
124      enddo
125!
126      return
127      end
128!
129!     subroutine df expressing the first order derivative as a function
130!     of time and the function itself
131!
132      subroutine df(x,u,uprime,rpar,nev)
133!
134      implicit none
135!
136      integer nev,i,j,id,iab,iaa,ibb,nev2p1,m
137!
138      real*8 rpar(*),x,u(*),uprime(*)
139!
140      iaa=4+nev*(nev+1)
141      ibb=iaa+nev
142      id=3+nev*nev
143!
144      do i=1,nev
145         uprime(i)=u(nev+i)
146         uprime(nev+i)=rpar(iaa+i)+x*rpar(ibb+i)
147     &             -rpar(id+i)*rpar(id+i)*u(i)
148     &             -(rpar(2)+rpar(3)*rpar(id+i)*rpar(id+i))*u(nev+i)
149!
150!        contribution of the dashpots
151!
152         do j=1,nev
153            uprime(nev+i)=uprime(nev+i)-rpar(3+(j-1)*nev+i)*u(nev+j)
154         enddo
155      enddo
156!
157      return
158      end
159!
160!     subroutine djac
161!
162      subroutine djac(x,u,pd,nrowpd,rpar,nev)
163!
164      implicit none
165!
166      integer nrowpd,nev,id,i,j
167!
168      real*8 rpar(*),x,u(*),pd(nrowpd,*)
169!
170      id=3+nev*nev
171!
172      do i=1,nev
173         pd(i,nev+i)=1.d0
174         pd(nev+i,i)=-rpar(id+i)*rpar(id+i)
175         pd(nev+i,nev+i)=-(rpar(2)+rpar(3)*rpar(id+i)*rpar(id+i))
176!
177!        contribution of the dashpots
178!
179         do j=1,nev
180            pd(nev+i,nev+j)=pd(nev+i,nev+j)-rpar(3+(j-1)*nev+i)
181         enddo
182      enddo
183!
184      return
185      end
186