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!     y=A*x for real sparse matrices (symmetric and non-symmetric)
20!
21!     storage of the matrix in a:
22!        - first the lower triangular terms
23!        - then, if the matrix is non-symmetric, the upper triangular terms
24!        - finally the diagonal terms
25!
26      subroutine matvec_struct(n,x,y,nelt,ia,ja,a,isym)
27      use omp_lib
28!
29      implicit none
30!
31      integer ia(*),ja(*),i,j,l,n,nelt,isym,nd,na
32      real*8 y(*),x(*),a(*)
33!
34!     number of off-diagonal terms
35!
36      nd=nelt-n
37!
38      if(isym.eq.0) then
39         na=nd/2
40!
41!        non-symmetric
42!
43!        diagonal terms
44!
45c$omp parallel default(none)
46c$omp& shared(n,x,a,y,nd,ja,ia,na)
47c$omp& private(i,l,j)
48c$omp do
49         do i=1,n
50            y(i)=a(nd+i)*x(i)
51         enddo
52c$omp end do
53!
54!        off-diagonal terms
55!
56!        number of upper triangular terms
57!
58c$omp do
59         do j=1,n
60            do l=ja(j),ja(j+1)-1
61               i=ia(l)
62c$omp atomic
63               y(i)=y(i)+a(l)*x(j)
64            enddo
65         enddo
66c$omp end do
67!
68c$omp do
69         do j=1,n
70            do l=ja(j),ja(j+1)-1
71               i=ia(l)
72               y(j)=y(j)+a(l+na)*x(i)
73            enddo
74         enddo
75c$omp end do
76c$omp end parallel
77!
78      else
79!
80!        symmetric
81!
82!        diagonal terms
83!
84c$omp parallel default(none)
85c$omp& shared(n,x,a,y,nd,ja,ia)
86c$omp& private(i,l,j)
87c$omp do
88         do i=1,n
89            y(i)=a(nd+i)*x(i)
90         enddo
91c$omp end do
92!
93!        off-diagonal terms
94!
95c$omp do
96         do j=1,n
97            do l=ja(j),ja(j+1)-1
98               i=ia(l)
99c$omp atomic
100               y(i)=y(i)+a(l)*x(j)
101            enddo
102         enddo
103c$omp end do
104!
105c$omp do
106         do j=1,n
107            do l=ja(j),ja(j+1)-1
108               i=ia(l)
109               y(j)=y(j)+a(l)*x(i)
110            enddo
111         enddo
112c$omp end do
113c$omp end parallel
114      endif
115!
116      return
117      end
118