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 frequencys(inpc,textpart,nmethod,
20     &  mei,fei,iperturb,istep,istat,n,iline,ipol,inl,
21     &  ipoinp,inp,ithermal,isolver,xboun,nboun,ipoinpc,
22     &  ipompc,labmpc,fmpc,ikmpc,ilmpc,nmpc,ier)
23!
24!     reading the input deck: *FREQUENCY
25!
26      implicit none
27!
28      logical global,cycmpcactive
29!
30      character*1 inpc(*)
31      character*20 solver,labmpc(*)
32      character*132 textpart(16)
33!
34      integer nmethod,mei(4),ncv,mxiter,istep,istat,iperturb(*),i,
35     &  nboun,ier,isolver,
36     &  n,key,iline,ipol,inl,ipoinp(2,*),inp(3,*),nev,ithermal(*),
37     &  ipoinpc(0:*),nmpcred,kflag,ipompc(*),ikmpc(*),ilmpc(*),nmpc
38!
39      real*8 fei(3),pi,fmin,fmax,tol,xboun(*),fmpc(*)
40!
41      global=.true.
42      cycmpcactive=.true.
43!
44      pi=4.d0*datan(1.d0)
45      mei(4)=0
46!
47!     defaults for fmin and fmax
48!
49      fmin=-1.d0
50      fmax=-1.d0
51!
52      if(istep.lt.1) then
53         write(*,*)
54     &      '*ERROR reading *FREQUENCY: *FREQUENCY can only be used'
55         write(*,*) '  within a STEP'
56         ier=1
57         return
58      endif
59!
60!     no heat transfer analysis
61!
62      if(ithermal(1).gt.1) then
63         ithermal(1)=1
64      endif
65!
66!     default solver
67!
68      solver='                    '
69      if(isolver.eq.0) then
70         solver(1:7)='SPOOLES'
71      elseif(isolver.eq.2) then
72         solver(1:16)='ITERATIVESCALING'
73      elseif(isolver.eq.3) then
74         solver(1:17)='ITERATIVECHOLESKY'
75      elseif(isolver.eq.4) then
76         solver(1:3)='SGI'
77      elseif(isolver.eq.5) then
78         solver(1:5)='TAUCS'
79      elseif(isolver.eq.7) then
80         solver(1:7)='PARDISO'
81      elseif(isolver.eq.8) then
82         solver(1:6)='PASTIX'
83      endif
84!
85      do i=2,n
86         if(textpart(i)(1:7).eq.'SOLVER=') then
87            read(textpart(i)(8:27),'(a20)') solver
88         elseif(textpart(i)(1:11).eq.'STORAGE=YES') then
89            mei(4)=1
90         elseif(textpart(i)(1:9).eq.'GLOBAL=NO') then
91            global=.false.
92         elseif(textpart(i)(1:15).eq.'CYCMPC=INACTIVE') then
93            cycmpcactive=.false.
94         else
95            write(*,*)
96     &        '*WARNING reading *FREQUENCY: parameter not recognized:'
97            write(*,*) '         ',
98     &                 textpart(i)(1:index(textpart(i),' ')-1)
99            call inputwarning(inpc,ipoinpc,iline,
100     &"*FREQUENCY%")
101         endif
102      enddo
103!
104      if(solver(1:7).eq.'SPOOLES') then
105         isolver=0
106      elseif(solver(1:16).eq.'ITERATIVESCALING') then
107         write(*,*) '*WARNING reading *FREQUENCY: the iterative scaling'
108         write(*,*) '         procedure is not available for frequency'
109         write(*,*) '         calculations; the default solver is used'
110      elseif(solver(1:17).eq.'ITERATIVECHOLESKY') then
111         write(*,*) '*WARNING reading *FREQUENCY: the iterative scaling'
112         write(*,*) '         procedure is not available for frequency'
113         write(*,*) '         calculations; the default solver is used'
114      elseif(solver(1:3).eq.'SGI') then
115         isolver=4
116      elseif(solver(1:5).eq.'TAUCS') then
117         isolver=5
118      elseif(solver(1:13).eq.'MATRIXSTORAGE') then
119         isolver=6
120      elseif(solver(1:7).eq.'PARDISO') then
121         isolver=7
122      elseif(solver(1:6).eq.'PASTIX') then
123         isolver=8
124      else
125         write(*,*) '*WARNING reading *FREQUENCY: unknown solver;'
126         write(*,*) '         the default solver is used'
127      endif
128!
129      if((isolver.eq.2).or.(isolver.eq.3)) then
130         write(*,*) '*ERROR reading *FREQUENCY: the default solver ',
131     & solver
132         write(*,*) '       cannot be used for frequency calculations '
133         ier=1
134         return
135      endif
136!
137      nmethod=2
138      if(iperturb(1).gt.1) iperturb(1)=0
139c      iperturb(2)=0
140!
141      if(isolver.ne.6) then
142         call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
143     &        ipoinp,inp,ipoinpc)
144         if((istat.lt.0).or.(key.eq.1)) then
145            write(*,*)
146     &       '*ERROR reading *FREQUENCY: definition not complete'
147            write(*,*) '  '
148            call inputerror(inpc,ipoinpc,iline,
149     &           "*FREQUENCY%",ier)
150            return
151         endif
152         read(textpart(1)(1:10),'(i10)',iostat=istat) nev
153         if(istat.gt.0) then
154            call inputerror(inpc,ipoinpc,iline,
155     &           "*FREQUENCY%",ier)
156            return
157         endif
158         if(nev.le.0) then
159            write(*,*) '*ERROR reading *FREQUENCY: less than 1 eigenvalu
160     &e requested'
161            ier=1
162            return
163         endif
164         tol=1.d-2
165         ncv=4*nev
166         ncv=ncv+nev
167         mxiter=1000
168         if(textpart(2)(1:1).ne.' ') then
169            read(textpart(2)(1:20),'(f20.0)',iostat=istat) fmin
170            if(istat.gt.0) then
171               call inputerror(inpc,ipoinpc,iline,
172     &              "*FREQUENCY%",ier)
173               return
174            endif
175         endif
176         if(textpart(3)(1:1).ne.' ') then
177            read(textpart(3)(1:20),'(f20.0)',iostat=istat) fmax
178            if(istat.gt.0) then
179               call inputerror(inpc,ipoinpc,iline,
180     &              "*FREQUENCY%",ier)
181               return
182            endif
183         endif
184!
185         mei(1)=nev
186         mei(2)=ncv
187         mei(3)=mxiter
188         fei(1)=tol
189         fei(2)=fmin
190         fei(3)=fmax
191!
192      else
193!
194!        matrix storage
195!
196         if(global) then
197            mei(1)=1
198         else
199            mei(1)=0
200         endif
201!
202         if(.not.cycmpcactive) then
203!
204!           remove the CYCLIC and SUBCYCLIC MPC's
205!
206            nmpcred=0
207!
208            kflag=2
209            call isortii(ilmpc,ikmpc,nmpc,kflag)
210!
211            do i=1,nmpc
212               if((labmpc(i)(1:6).eq.'CYCLIC').or.
213     &            (labmpc(i)(1:9).eq.'SUBCYCLIC')) cycle
214               nmpcred=nmpcred+1
215               ipompc(nmpcred)=ipompc(i)
216               labmpc(nmpcred)=labmpc(i)
217               fmpc(nmpcred)=fmpc(i)
218               ikmpc(nmpcred)=ikmpc(i)
219               ilmpc(nmpcred)=nmpcred
220            enddo
221!
222            nmpc=nmpcred
223            call isortii(ikmpc,ilmpc,nmpc,kflag)
224         endif
225      endif
226!
227!     removing nonzero boundary conditions
228!
229      do i=1,nboun
230         xboun(i)=0.d0
231      enddo
232!
233      call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
234     &     ipoinp,inp,ipoinpc)
235!
236      return
237      end
238
239