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!      INPUT Parameters:
20!      neq:      number of equations, length of eigenvectors
21!                (non-cyclic-symmetry)
22!      nev:      number of eigenvectors
23!      neqact:   length of eigenvectors (cyclic symmetry)
24!      z:        Matrix of eigenvectors out of frequency analysis
25!                storage: for each eigenvalue: first complete real part,
26!                         then complete imaginary part of eigenvector
27!                         (imaginary part only for cyclic symmetry)
28!      zz:       Matrix of eigenvectors out of complex frequency analysis
29!                storage: for each eigenvalue: first complete real part,
30!                         then complete imaginary part of eigenvector
31!      istartnmd:first number of nev of nodal diameter l
32!      iendnmd:  last number ov nev of nodal diameter l
33!      nmd:      Number of nodal diameters (cyclic symmetry)
34!
35!      OUTPUT Parameters
36!      xmac:     Matrix that contains absolute MAC-values of all vectors
37!                combined
38!                with all vectors; in case of cyclic symmetry: only eigenvectors
39!                of the same nodal diameter are taken in account
40!      xmaccpx:  Matrix that contains the complex MAC-values of all vectors
41!                combined with all vectors; in case of cyclic symmetry:
42!                only eigenvectors of the same nodal diameter are taken
43!                in account
44!
45      subroutine calcmac(neq,z,zz,nev,xmac,xmaccpx,istartnmd,
46     &     iendnmd,nmd,cyclicsymmetry,neqact,bett,betm,nevcomplex)
47!
48!     calculates the Modal Assurance Criterium MAC=<z,zz>/(||z||*||zz||)
49!
50      implicit none
51!
52      integer neq,nev,i,j,k,l,istartnmd(*),iendnmd(*),nmd,
53     &     cyclicsymmetry,neqact,nevcomplex
54!
55      real*8 bett(*),betm(*),xmac(nev,*),z(neq,*),zz(2*neqact,*)
56!
57      complex*16 xmaccpx(nev,*)
58!
59      if(cyclicsymmetry.eq.0)then
60!     size of vectors
61         do i=1,nev
62            do k=1,neq
63               bett(i)=bett(i)+z(k,i)**2
64            enddo
65         enddo
66         do i=1,nevcomplex
67            do k=1,neq
68               betm(i)=betm(i)+zz(k,i)**2+zz(k+neq,i)**2
69            enddo
70         enddo
71!
72         do i=1,nev
73            bett(i)=dsqrt(bett(i))
74         enddo
75         do i=1,nevcomplex
76            betm(i)=dsqrt(betm(i))
77         enddo
78!     Calculation of MAC
79         do i=1,nev
80            do j=1,nevcomplex
81               do k=1,neq
82                  xmaccpx(i,j)=xmaccpx(i,j)+z(k,i)
83     &                 *(zz(k,j)+zz(k+neq,j)*(0.d0,1.d0))
84               enddo
85               xmac(i,j)=cdabs(xmaccpx(i,j))
86               xmac(i,j)=xmac(i,j)/bett(i)/betm(j)
87            enddo
88         enddo
89!
90!     Cyclic Symmetry
91!     size of vectors
92!
93      else
94         do i=1,nev
95            do k=1,neqact
96               bett(i)=bett(i)+z(k,i)**2+z(k+neqact,i)**2
97               betm(i)=betm(i)+zz(k,i)**2+zz(k+neqact,i)**2
98            enddo
99            bett(i)=dsqrt(bett(i))
100            betm(i)=dsqrt(betm(i))
101         enddo
102!     Calculation of MAC
103         do l=1,nmd
104            do i=istartnmd(l),iendnmd(l)
105               do j=istartnmd(l),iendnmd(l)
106                  xmac(i,j)=0
107                  do k=1,neqact,2
108                     xmaccpx(i,j)=xmaccpx(i,j)+
109     &                    (z(k,i)-z(k+neqact,i)*(0.d0,1.d0))*
110     &                    (zz(k,j)+zz(k+neqact,j)*(0.d0,1.d0))
111                  enddo
112                  xmac(i,j)=cdabs(xmaccpx(i,j))
113                  xmac(i,j)=(xmac(i,j))/bett(i)/betm(j)
114               enddo
115            enddo
116         enddo
117      endif
118!
119      return
120      end
121
122