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