DOUG 0.2

SpMtx_op_Ax.f90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 !!----------------------------------------------------------
00023 !!Operation Ax for Sparse Matrix
00024 !! A - sparse Matrix
00025 !! x - vector
00026 !!Include:
00027 !!        simple operation Ax (for not arranged matrix)
00028 !!        operation Ax for arranged matrix (fast)
00029 !!----------------------------------------------------------
00030 Module SpMtx_op_Ax
00031   use RealKind
00032   use SpMtx_class
00033   use SpMtx_arrangement
00034   Implicit None
00035 CONTAINS
00036 !----------------------------------------------------------
00037 !Simple Operation Ax
00038 !    Arguments:
00039 !           A - sparse matrix (type SpMtx)
00040 !           x - vector
00041 !      dozero - .TRUE. : y=0
00042 !    Result: Vector y=Ax
00043 !----------------------------------------------------------
00044   subroutine SpMtx_Ax(y,A,x,dozero,transp)
00045     Implicit None
00046     type(SpMtx), intent(in)                :: A        !sparse matrix
00047     real(kind=rk),dimension(:),pointer     :: x        !vector (in)
00048     real(kind=rk),dimension(:),pointer     :: y        !result vector
00049     integer                                :: i,ii,j,arrtype,ncols,nrows
00050     logical, optional, intent(in)          :: dozero   !
00051     logical, optional, intent(in)          :: transp   !
00052     integer,dimension(:),pointer           :: indi,indj
00053     !- - - - - - - - - - - - - - - - - - - - - - - - - - -
00054     if (present(dozero)) then
00055       if (dozero) y=0.0_rk
00056     endif
00057     if (present(transp).and.(transp)) then
00058       indi=>A%indj
00059       indj=>A%indi
00060       if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then
00061         arrtype=D_SpMtx_ARRNG_COLS
00062       elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then
00063         arrtype=D_SpMtx_ARRNG_ROWS
00064       else
00065         arrtype=A%arrange_type
00066       endif
00067       nrows=A%ncols
00068       ncols=A%nrows
00069     else
00070       indi=>A%indi
00071       indj=>A%indj
00072       arrtype=A%arrange_type
00073       nrows=A%nrows
00074       ncols=A%ncols
00075     endif
00076 
00077     if (arrtype==D_SpMtx_ARRNG_NO) then
00078       do j=1,A%nnz
00079         i=indi(j)
00080         y(i)=y(i)+A%val(j)*x(indj(j))
00081       enddo
00082     elseif (arrtype==D_SpMtx_ARRNG_ROWS) then
00083       do i=1,nrows
00084         do j=A%M_bound(i),A%M_bound(i+1)-1
00085           y(i)=y(i)+A%val(j)*x(indj(j))
00086         end do
00087       end do
00088     elseif (arrtype==D_SpMtx_ARRNG_COLS) then
00089 !rite(stream,*)'ncols:',ncols
00090 !rite(stream,*)'size(y):',size(y),y(1:5)
00091 !rite(stream,*)'y(1:5):',y(1:5)
00092       do j=1,ncols
00093         do i=A%M_bound(j),A%M_bound(j+1)-1
00094           ii=indi(i)
00095           y(ii)=y(ii)+A%val(i)*x(j)
00096         end do
00097       end do
00098     else
00099       call DOUG_abort('[SpMtx_Ax] : Incorrect matrix arrange type.', -1)
00100    endif
00101   end subroutine SpMtx_Ax
00102 !----------------------------------------------------------
00103 !Operation Ax for arranged matrix
00104 !  If Matrix A does not be arranged...program does that.
00105 !    Arguments:
00106 !           A - sparse matrix (arranged by rows)
00107 !           x - vector
00108 !    Result: Vector Sx=Ax
00109 !----------------------------------------------------------
00110   Function SpMtx_arrangedAx(A,x) result(Sx)
00111     Implicit None
00112     type(SpMtx), intent(in out)        :: A    !sparse matrix (in)
00113     real(kind=rk), intent(in), dimension(:):: x    !vector (in)
00114     real(kind=rk), dimension(:), pointer   :: Sx   !result vector
00115     integer                                :: i, j !counters
00116     !- - - - - - - - - - - - - - - - - - - - - - - - -
00117     if (A%arrange_type /= 1) call SpMtx_arrange(A) !arrange if nessesary
00118     if (size(x) /= A%ncols) print*, "ERROR: Ax dim(A) /= dim(x)"
00119     allocate(Sx(A%nrows)); Sx=0.
00120     do i=1,A%nrows                            !calculate result
00121       do j=A%M_bound(i),A%M_bound(i+1)-1
00122         Sx(i)=Sx(i)+A%val(j)*x(A%indj(j))
00123       end do
00124     end do
00125   End Function SpMtx_arrangedAx
00126 End Module SpMtx_op_Ax
00127 !----------------------------------------------------------
00128 !$Log: SpMtx_op_Ax.f90,v $
00129 !Revision 1.6  2004/04/30 08:57:24  elmo
00130 !Formatting improved
00131 !
00132 !Revision 1.5  2004/03/08 07:37:15  elmo
00133 !Added files for AMG and some test files.
00134 !
00135 !Revision 1.4  2003/12/04 08:08:25  eero
00136 !Made some changes to get PCG (version 2) running correctly with ILU on SUN.
00137 !
00138 !Revision 1.3  2003/11/03 08:43:39  elmo
00139 !Changed operation SpMtx_Ax, added files for block_m type
00140 !
00141 !Revision 1.2  2003/10/30 14:31:54  elmo
00142 !Fixed syntax bug
00143 !
00144 !Revision 1.1  2003/10/30 14:21:57  elmo
00145 !Added files SpMtx_arrange.f90 SpMtx_Ax.f90
00146 !
00147 !----------------------------------------------------------