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 forcesolve(zc,nev,a,b,x,eiga,eigb,eigxx,iter,d,
20     &  neq,z,istartnmd,iendnmd,nmd,cyclicsymmetry,neqact,
21     &  igeneralizedforce)
22!
23!     solves for the complex eigenfrequencies due to Coriolis
24!     forces
25!
26      implicit none
27!
28      logical wantx
29!
30      integer nev,neq,iter(*),i,j,k,l,istartnmd(*),iendnmd(*),nmd,
31     &  cyclicsymmetry,neqact,igeneralizedforce
32!
33      real*8 z(neq,*),d(*)
34!
35      complex*16 a(nev,*),b(nev,*),x(nev,*),eiga(*),eigb(*),eigxx(*),
36     &  zc(neqact,*)
37!
38      if(igeneralizedforce.eq.0) then
39!
40!        no generalized force: multiplication with the eigenmodes
41!        is necessary
42!
43         if(cyclicsymmetry.eq.0) then
44            do i=1,nev
45               do j=1,nev
46                  do k=1,neq
47                     a(i,j)=a(i,j)+z(k,i)*zc(k,j)
48                  enddo
49               enddo
50               write(*,*)
51     &              'aerodynamic stiffness/structural stiffness = ',
52     &              a(i,i)/d(i)
53               a(i,i)=a(i,i)+d(i)
54               b(i,i)=(1.d0,0.d0)
55            enddo
56         else
57!
58!     cyclic symmetry
59!
60            do l=1,nmd
61               do i=istartnmd(l),iendnmd(l)
62                  do j=istartnmd(l),iendnmd(l)
63                     do k=1,neqact
64                        a(i,j)=a(i,j)+z(k,i)*zc(k,j)-
65     &                       z(k+neqact,i)*zc(k,j)*(0.d0,1.d0)
66                     enddo
67                  enddo
68                  write(*,*)
69     &                 'aerodynamic stiffness/structural stiffness = ',
70     &                 a(i,i)/d(i)
71                  a(i,i)=a(i,i)+d(i)
72                  b(i,i)=(1.d0,0.d0)
73               enddo
74            enddo
75         endif
76      else
77!
78!        generalized force: the a-matrix is (apart from the diagonal)
79!        known
80!
81         if(cyclicsymmetry.eq.0) then
82            do i=1,nev
83               write(*,*)
84     &              'aerodynamic stiffness/structural stiffness = ',
85     &              a(i,i)/d(i)
86               a(i,i)=a(i,i)+d(i)
87               b(i,i)=(1.d0,0.d0)
88            enddo
89         else
90!
91!     cyclic symmetry
92!
93            do l=1,nmd
94               do i=istartnmd(l),iendnmd(l)
95                  write(*,*)
96     &                 'aerodynamic stiffness/structural stiffness = ',
97     &                 a(i,i)/d(i)
98                  a(i,i)=a(i,i)+d(i)
99                  b(i,i)=(1.d0,0.d0)
100               enddo
101            enddo
102         endif
103      endif
104!
105      wantx=.true.
106!
107!     solving for the complex eigenvalues
108!
109      call dlzhes(nev,a,nev,b,nev,x,nev,wantx)
110      call dlzit(nev,a,nev,b,nev,x,nev,wantx,iter,eiga,eigb)
111!
112      do i=1,nev
113         if(iter(i).eq.-1) then
114            write(*,*) '*ERROR in forcesolve: fatal error'
115            write(*,*) '       in dlzit'
116            call exit(201)
117         elseif(cdabs(eigb(i)).lt.1.d-10) then
118            write(*,*) '*ERROR in forcesolve: eigenvalue'
119            write(*,*) '       out of bounds'
120            call exit(201)
121         else
122            eigxx(i)=cdsqrt(eiga(i)/eigb(i))
123         endif
124      enddo
125!
126      return
127      end
128
129
130