fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

New files from Thomas.

Browse Source
tamestoy
Thomas Amestoy 11 years ago
parent d864f5b60e
commit 9068df9228
11 changed files with 2342 additions and 2 deletions
Show Stats Download Patch File Download Diff File

11
test/kernel/Makefile
Unescape Escape View File

@ -10,12 +10,11 @@ LDLIBS=$(PSBLDLIBS)
FINCLUDES=$(FMFLAG)$(INCDIR) $(FMFLAG).
DTOBJS=d_file_spmv.o
STOBJS=s_file_spmv.o
DPGOBJS=pdgenspmv.o
PWFOBJS = power_file.o
PWFROBJS = power_file_real.o
PWFROBJS = power_file_real.o #eig_mod.o
ARNOBJS =arnoldi_file.o
ARNROBJS =arnoldi_file_real.o
SIROBJS=shift_invert_real.o
@ -60,6 +59,10 @@ analyse: analyse.o
$(F90LINK) $(LOPT) analyse.o -o analyse $(PSBLAS_LIB) $(LDLIBS)
/bin/mv analyse /server/tamestoy/TEMP/MaxEigenAdj/
an_lapl: an_lapl.o
$(F90LINK) $(LOPT) an_lapl.o -o an_lapl $(PSBLAS_LIB) $(LDLIBS)
/bin/mv an_lapl /server/tamestoy/TEMP/SCALING-BA
laplacian: laplacian.o
$(F90LINK) $(LOPT) laplacian.o -o laplacian $(PSBLAS_LIB) $(LDLIBS)
/bin/mv laplacian $(EXEDIR)
@ -67,6 +70,10 @@ laplacian: laplacian.o
test_chseqr: test_chseqr.o
$(F90LINK) $(LOPT) test_chseqr.o -o test_chseqr $(PSBLAS_LIB) $(LDLIBS)
#eig_mod.o: eig_mod.f90
#$(F90LINK) $(LOPT) -c eig_mod.f90 $(PSBLAS_LIB) $(LDLIBS)
#gfortran -c eig_mod.f90
clean:
/bin/rm -f *.o
#/bin/rm -f $(DBOBJSS) $(DBOBJS) $(PWOBJS) $(PWFOBJS) $(DTOBJS) $(STOBJS) $(ARNOBJS)

38
test/kernel/adj_read.f90
Unescape Escape View File

@ -0,0 +1,38 @@
subroutine adj_read (a,filename,desc_a,info)
type(psb_dspmat_type), intent (out) :: a
character(len=20) :: filename
type(psb_desc_type):: desc_a
integer (psb_ipk_) :: info
implicit none
integer :: i,j,nnzero,nbCol
integer :: unitFile, iError,line
integer(psb_ipk_), allocatable :: fileContent(:,:),ia(:),ja(:)
real (psb_dpk_), allocatable :: val(:)
integer(psb_ipk_), allocatable :: size_mat(:)
nbCol = 2
allocate (size_mat(nbCol))
unitFile = 1
open(UNIT=unitFile, FILE=filename, FORM="FORMATTED", STATUS="OLD",
ACTION="READ")
nnzero = size_mat(2)
allocate (fileContent(nnzeros,nbCol))
do line = 1,nnzero
read(unitFile, *) fileContent(line,1:nbCol)
end do saveNodes
close(UNIT=unitFile)
allocate(ia(nnzero),ja(nnzero),val(nnzero))
do i=1,nnzero
ia(i)=fileContent(i,1)
ja(i)=filecontent(i,2)
val(i)=1.0
end do
call psb_spins(nnzero, ia, ja, val, a, desc_a, info)
end subroutine adj_read

82
test/kernel/an_lapl.f90
Unescape Escape View File

@ -0,0 +1,82 @@
program analyse
use psb_base_mod
use psb_util_mod
implicit none
character(len=30) :: lapl_file,max_file,res_file,analyse_file,mat_file
integer (psb_ipk_) :: i,nb_mat, n
real(psb_dpk_) ::tab(18),eig1,eig2,eig0,eigmin,eigmax,delta,delta2,&
&deltamax,timemin,timemax,t1,t2,time0,t0
read(psb_inp_unit,*) lapl_file
read(psb_inp_unit,*) max_file
read(psb_inp_unit,*) res_file
read(psb_inp_unit,*) analyse_file
read(psb_inp_unit,*) nb_mat
open(15,FILE=lapl_file,action="read")
open(13,FILE=max_file,action="read")
open(14,file=res_file,action="read")
do i=1,21
read(15,*)
end do
do i=1,22
read(13,*)
end do
open (16, file=analyse_file,position="append",action="write")
n=0
delta=0
delta2=0
deltamax=0
timemax=0
timemin=0
time0=0
do i=1,nb_mat
read(15,*)tab(1:18)
read(14,*)mat_file,eig1,eig2,t1
read(14,*)eigmin,eig0,t2
if(abs(tab(6)-eigmin)>0.001)then
n=n+1
write(16,'("smallest of : ",F20.2,F20.2)')tab(2),tab(3)
write(16,'(F20.7,F20.7)')tab(6),eigmin
end if
if(abs(tab(16)-eig2)>0.01)then
n=n+1
write(16,'("lambda N-1 of : ",F20.2,F20.2)')tab(2),tab(3)
write(16,'(F20.7,F20.7)')tab(16),eig2
end if
if(abs(tab(18)-eig1)>0.01)then
n=n+1
write(16,'("biggest of : ", F20.2,F20.2)')tab(2),tab(3)
write(16,'(F20.7,F20.7)')tab(18),eig1
end if
timemax=timemax+t1
timemin=timemin+t2
delta2=delta2+abs(1-eigmin/tab(6))
delta=delta+abs(1-eig2/tab(16))+abs(1-eig1/tab(18))
read(13,*)tab(1:7)
read(14,*)mat_file,eigmax,t0
if(abs(tab(4)-eigmax)>0.01)then
write(16,'(F20.2,F20.2)')tab(2),tab(3)
write(16,'(F20.6,F20.6)')tab(4),eigmax
endif
time0=time0+t0
deltamax=deltamax+abs(1-eigmax/tab(4))
end do
write(16,'("errors on ",i20," eigenvalues")')n
write(16,'("gap average for lapl max",g20.5)')delta/(2*nb_mat)
write(16,'("gap average for lapl min",g20.5)')delta2/(nb_mat)
write(16,'("gap average for max",g20.5)')deltamax/(nb_mat)
write(16,'("time average lapl max",g20.5)')timemax/(nb_mat)
write(16,'("time average lapl min",g20.5)')timemin/(nb_mat)
write(16,'("time average A max",g20.5)')time0/(nb_mat)
close(15)
close(14)
close(13)
close(16)
end program analyse

43
test/kernel/analyse.f90
Unescape Escape View File

@ -0,0 +1,43 @@
program analyse
use psb_base_mod
use psb_util_mod
implicit none
character(len=30) :: given_file,res_file,analyse_file,mat_file
integer (psb_ipk_) :: i,nb_mat, n
real(psb_dpk_) :: tab(7),eig
read(psb_inp_unit,*) given_file
read(psb_inp_unit,*) res_file
read(psb_inp_unit,*) analyse_file
read(psb_inp_unit,*) nb_mat
open(15,FILE=given_file,action="read")
open(14,file=res_file,action="read")
do i=1,22
read(15,*)
end do
open (16, file=analyse_file,position="append",action="write")
n=0
do i=1,nb_mat
read(15,*)tab(1:7)
read(14,*)mat_file,eig
if(tab(4)-eig==0)then
write(psb_out_unit,'("on a gagne !")')
n=n+1
else
write(psb_out_unit,'("perdu")')
write(16,'(F20.2,F20.2)')tab(2),tab(3)
write(16,'(F20.7,F20.7)')tab(4),eig
end if
end do
write(psb_out_unit,'(i20)')n
close(14)
close(15)
close(16)
end program analyse

343
test/kernel/arnoldi_file.f90
Unescape Escape View File

@ -0,0 +1,343 @@
!!$
!!$ Parallel Sparse BLAS version 3.1
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
program arnoldi_file
use psb_base_mod
use psb_util_mod
implicit none
! input parameters
character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
! sparse matrices
type(psb_zspmat_type) :: a, aux_a
! dense matrices
complex(psb_dpk_), allocatable, target :: aux_b(:,:), d(:)
complex(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:)
complex(psb_dpk_), allocatable, target :: H(:,:),eig(:),work(:),Z(:,:)
integer, allocatable :: indexes(:)
type(psb_z_vect_type), allocatable, target :: V(:)
complex(psb_dpk_), pointer :: b_col_glob(:)
type(psb_z_vect_type) :: b_col, x_col, r_col
! communications data structure
type(psb_desc_type):: desc_a
integer(psb_ipk_) :: ictxt, iam, np
! solver paramters
integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
& methd, istopc, irst, nr
integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes
real(psb_dpk_) :: err, eps,cond
character(len=5) :: afmt
character(len=20) :: name
character(len=2) :: filefmt
integer(psb_ipk_), parameter :: iunit=12
integer(psb_ipk_) :: times=0
integer(psb_ipk_) :: iparm(20)
! other variables
integer(psb_ipk_) :: i,info,j,m_problem
integer(psb_ipk_) :: internal, m,ii,nnzero, dim_H, alloc_stat
real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
&scale,resmx,resmxp, flops, bdwdth
real(psb_dpk_) :: tt1, tt2, tflops
real (psb_dpk_) :: norm
complex (psb_dpk_) :: dotprod
integer(psb_ipk_) :: nrhs, nrow, n_row, nv, ne
integer(psb_ipk_), allocatable :: ivg(:), ipv(:)
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ictxt)
stop
endif
name='arnoldi_file'
if(psb_get_errstatus() /= 0) goto 9999
info=psb_success_
call psb_set_errverbosity(2)
!
! Hello world
!
if (iam == psb_root_) then
write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
write(*,*) 'This is the ',trim(name),' sample program'
read(psb_inp_unit,*) mtrx_file
read(psb_inp_unit,*) filefmt
read(psb_inp_unit,*) ipart
read(psb_inp_unit,*) dim_H
write (psb_out_unit, '("Hessenberg matrix dim is ",i20)') dim_H
end if
call psb_bcast(ictxt,mtrx_file)
call psb_bcast(ictxt,filefmt)
call psb_bcast(ictxt,ipart)
call psb_bcast(ictxt,dim_H)
rhs_file = 'NONE'
afmt = 'CSR'
call psb_barrier(ictxt)
t1 = psb_wtime()
! read the input matrix to be processed and (possibly) the rhs
nrhs = 1
if (iam==psb_root_) then
select case(psb_toupper(filefmt))
case('MM')
! For Matrix Market we have an input file for the matrix
! and an (optional) second file for the RHS.
call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
if (info == psb_success_) then
if (rhs_file /= 'NONE') then
call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
end if
end if
case ('HB')
! For Harwell-Boeig we have a single file which may or may not
! contain an RHS.
call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
case default
info = -1
write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
end select
if (info /= psb_success_) then
write(psb_err_unit,*) 'Error while reading input matrix '
call psb_abort(ictxt)
end if
m_problem = aux_a%get_nrows()
call psb_bcast(ictxt,m_problem)
! At this point aux_b may still be unallocated
if (psb_size(aux_b,dim=1)==m_problem) then
! if any rhs were present, broadcast the first one
write(psb_err_unit,'("Ok, got an rhs ")')
b_col_glob =>aux_b(:,1)
else
write(psb_out_unit,'("Generating an rhs...")')
write(psb_out_unit,'(" ")')
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob => aux_b(:,1)
do i=1, m_problem
b_col_glob(i) = 1.d0
enddo
endif
call psb_bcast(ictxt,b_col_glob(1:m_problem))
else
call psb_bcast(ictxt,m_problem)
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob =>aux_b(:,1)
call psb_bcast(ictxt,b_col_glob(1:m_problem))
end if
! switch over different partition types
if (ipart == 0) then
call psb_barrier(ictxt)
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
allocate(ivg(m_problem),ipv(np))
do i=1,m_problem
call part_block(i,m_problem,np,ipv,nv)
ivg(i) = ipv(1)
enddo
call psb_matdist(aux_a, a,ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else if (ipart == 2) then
if (iam==psb_root_) then
write(psb_out_unit,'("Partition type: graph")')
write(psb_out_unit,'(" ")')
! write(psb_err_unit,'("Build type: graph")')
call build_mtpart(aux_a,np)
endif
call psb_barrier(ictxt)
call distr_mtpart(psb_root_,ictxt)
call getv_mtpart(ivg)
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block)
end if
allocate(H(dim_H,dim_H),stat = alloc_stat)
do i=1,dim_h
do j=1,dim_H
H(i,j)=zzero
enddo
enddo
allocate(V(dim_H+1),stat = alloc_stat)
do i=1,dim_H+1
call psb_geall(V(i),desc_a,info)
call psb_geasb(V(i),desc_a,info)
enddo
call V(1)%set(zone)
t2 = psb_wtime() - t1
call psb_amx(ictxt, t2)
if (iam==psb_root_) then
write(psb_out_unit,'(" ")')
write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
write(psb_out_unit,'(" ")')
end if
call psb_barrier(ictxt)
t1 = psb_wtime()
norm = psb_genrm2(V(1),desc_a,info)
if (iam==psb_root_) write(psb_out_unit,'("norma iniziale :"F20.3)')norm
! H(2,1)=dcmplx(norm,0.0)
H(2,1)=norm*zone
norm = 1/norm
!normalisation of V(1)
call psb_geaxpby(zzero,V(1),norm*zone,V(1),desc_a, info)
do i=2,dim_H+1
call psb_spmm(zone,a,V(i-1),zzero,V(i),desc_a,info,'n') !we do V(i)=a*V(i-1)
! Gram-Schmitt's reorthogonalisation
do j=1,i-1
dotprod= psb_gedot(V(i),V(j),desc_a,info) ! dotprod = (V(i) dot V(j))
call psb_geaxpby(-dotprod,V(j),zone,V(i),desc_a, info) !V(i)=V(i)-V(j)*dotprod
if(iam==psb_root_) write(*,'("dotprod : "i5, g20.4,g20.4)') i,real(dotprod),aimag(dotprod)
H(j,i-1)=dotprod
end do
norm = psb_genrm2(V(i),desc_a,info)
if (iam==psb_root_) then
write(psb_out_unit,'("norma finale :"i20,F20.3)')i,norm
write(psb_out_unit,'("")')
end if
if (i .ne. dim_H+1) then
H(i,i-1)=cmplx(norm,0.0)
!H(i,i-1)=norm*zone
endif
norm=1/norm
call psb_geaxpby(zzero,V(i),norm*zone,V(i),desc_a, info)
enddo
do i=2,dim_H
if (iam==psb_root_) write(*,'("basse diagonale de H : "i5, g20.4,g20.4)') i,real(H(i,i-1)),aimag(H(i,i-1))
enddo
write(psb_out_unit,'("")')
allocate(eig(dim_H),work(dim_h),stat = info)
allocate(Z(dim_H,dim_H),stat = alloc_stat)
call ZHSEQR('E','N',dim_H,1,dim_H,H,dim_H,eig,Z,dim_H,work,3*dim_H,info)
!sort H's eigenvalues
allocate(indexes(1:dim_H))
call psb_qsort(eig,indexes,psb_alsort_up_,psb_sort_ovw_idx_)
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
call psb_amx(ictxt,t2)
nr = desc_a%get_global_rows()
annz = a%get_nzeros()
amatsize = psb_sizeof(a)
descsize = psb_sizeof(desc_a)
call psb_sum(ictxt,annz)
call psb_sum(ictxt,amatsize)
call psb_sum(ictxt,descsize)
if (iam==psb_root_) then
flops = 2.d0*times*annz
tflops=flops
write(psb_out_unit,'("Matrix: ",a)') mtrx_file
write(psb_out_unit,'("Test on : ",i20," processors")') np
write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
write(*,'("valeurs propres de H : ")')
OPEN(unit=10, file=mtrx_file(1:3)//'.txt')
do i=dim_H/3,dim_H
write(psb_out_unit,'(g20.4,g20.4)')real(eig(i)),aimag(eig(i))
write(10,'(g20.4,g20.4)')real(eig(i)),aimag(eig(i))
enddo
CLOSE(10)
tflops = tflops / (tt2)
!
! This computation is valid for CSR
!
nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_int)+ &
& annz*(psb_sizeof_sp + psb_sizeof_int)
bdwdth = times*nbytes/(t2*1.d6)
bdwdth = times*nbytes/(tt2*1.d6)
end if
call psb_gefree(b_col, desc_a,info)
call psb_gefree(x_col, desc_a,info)
call psb_spfree(a, desc_a,info)
call psb_cdfree(desc_a,info)
do i=1,dim_H
call psb_gefree(V(i),desc_a,info)
end do
DEALLOCATE (H)
DEALLOCATE (eig,V)
DEALLOCATE (work)
DEALLOCATE (Z)
9999 continue
if(info /= 0) then
call psb_error(ictxt)
end if
call psb_exit(ictxt)
stop
end program arnoldi_file

409
test/kernel/arnoldi_file_real.f90
Unescape Escape View File

@ -0,0 +1,409 @@
!!$
!!$ Parallel Sparse BLAS version 3.1
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
program arnoldi_file
use psb_base_mod
use psb_util_mod
implicit none
! input parameters
character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
! sparse matrices
type(psb_dspmat_type) :: a, b, aux_a
! dense matrices
real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:)
real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:)
complex(psb_dpk_), allocatable, target :: H(:,:),eig(:),work(:),Z(:,:)
integer, allocatable :: indexes(:)
type(psb_d_vect_type), allocatable, target :: V(:)
real(psb_dpk_), pointer :: b_col_glob(:)
type(psb_d_vect_type) :: b_col, x_col, r_col
! communications data structure
type(psb_desc_type):: desc_a
integer(psb_ipk_) :: ictxt, iam, np
! solver paramters
integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
& methd, istopc, irst, nr
integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes
real(psb_dpk_) :: err, eps,cond
character(len=5) :: afmt
character(len=20) :: name
character(len=2) :: filefmt
integer(psb_ipk_), parameter :: iunit=12
integer(psb_ipk_) :: times=0
integer(psb_ipk_) :: iparm(20)
! other variables
integer(psb_ipk_) :: i,info,j,m_problem
integer(psb_ipk_) :: internal, m,ii,nnzero, dim_H, alloc_stat
real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
&scale,resmx,resmxp, flops, bdwdth
real(psb_dpk_) :: tt1, tt2, tflops
real (psb_dpk_) :: norm
real (psb_dpk_) :: dotprod
integer(psb_ipk_) :: nrhs, nrow, n_row, nv, ne
integer(psb_ipk_), allocatable :: ivg(:), ipv(:)
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ictxt)
stop
endif
name='arnoldi_file'
if(psb_get_errstatus() /= 0) goto 9999
info=psb_success_
call psb_set_errverbosity(2)
!
! Hello world
!
if (iam == psb_root_) then
!write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
!write(*,*) 'This is the ',trim(name),' sample program'
read(psb_inp_unit,*) mtrx_file
read(psb_inp_unit,*) filefmt
read(psb_inp_unit,*) ipart
read(psb_inp_unit,*) dim_H
end if
call psb_bcast(ictxt,mtrx_file)
call psb_bcast(ictxt,filefmt)
call psb_bcast(ictxt,ipart)
call psb_bcast(ictxt,dim_H)
rhs_file = 'NONE'
afmt = 'CSR'
call psb_barrier(ictxt)
t1 = psb_wtime()
! read the input matrix to be processed and (possibly) the rhs
nrhs = 1
if (iam==psb_root_) then
select case(psb_toupper(filefmt))
case('MM')
! For Matrix Market we have an input file for the matrix
! and an (optional) second file for the RHS.
call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
if (info == psb_success_) then
if (rhs_file /= 'NONE') then
call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
end if
end if
case ('HB')
! For Harwell-Boeig we have a single file which may or may not
! contain an RHS.
call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
case ('AD')
call adj_read(aux_a,mtrx_file,iunit,desc_a,info)
case default
info = -1
write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
end select
if (info /= psb_success_) then
write(psb_err_unit,*) 'Error while reading input matrix '
call psb_abort(ictxt)
end if
m_problem = aux_a%get_nrows()
call psb_bcast(ictxt,m_problem)
! At this point aux_b may still be unallocated
if (psb_size(aux_b,dim=1)==m_problem) then
! if any rhs were present, broadcast the first one
write(psb_err_unit,'("Ok, got an rhs ")')
b_col_glob =>aux_b(:,1)
else
! write(psb_out_unit,'("Generating an rhs...")')
! write(psb_out_unit,'(" ")')
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob => aux_b(:,1)
do i=1, m_problem
b_col_glob(i) = 1.d0
enddo
endif
call psb_bcast(ictxt,b_col_glob(1:m_problem))
else
call psb_bcast(ictxt,m_problem)
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob =>aux_b(:,1)
call psb_bcast(ictxt,b_col_glob(1:m_problem))
end if
! switch over different partition types
if (ipart == 0) then
call psb_barrier(ictxt)
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
allocate(ivg(m_problem),ipv(np))
do i=1,m_problem
call part_block(i,m_problem,np,ipv,nv)
ivg(i) = ipv(1)
enddo
call psb_matdist(aux_a, a,ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else if (ipart == 2) then
if (iam==psb_root_) then
!write(psb_out_unit,'("Partition type: graph")')
!write(psb_out_unit,'(" ")')
! write(psb_err_unit,'("Build type: graph")')
call build_mtpart(aux_a,np)
endif
call psb_barrier(ictxt)
call distr_mtpart(psb_root_,ictxt)
call getv_mtpart(ivg)
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else
!if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block)
end if
call lapl(a,b)
allocate(H(dim_H,dim_H),stat = alloc_stat)
do i=1,dim_h
do j=1,dim_H
H(i,j)=zzero
enddo
enddo
allocate(V(dim_H+1),stat = alloc_stat)
do i=1,dim_H+1
call psb_geall(V(i),desc_a,info)
call psb_geasb(V(i),desc_a,info)
enddo
call V(1)%set(done)
t2 = psb_wtime() - t1
call psb_amx(ictxt, t2)
! if (iam==psb_root_) then
!write(psb_out_unit,'(" ")')
!write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
!write(psb_out_unit,'(" ")')
!end if
call psb_barrier(ictxt)
norm = psb_norm2(V(1),desc_a,info)
H(2,1)=cmplx(norm,0.0)
norm = 1/norm
!normalisation of V(1)
call psb_geaxpby(dzero,V(1),norm,V(1),desc_a, info)
do i=2,dim_H+1
call psb_spmm(done,b,V(i-1),dzero,V(i),desc_a,info,'n') !we do V(i)=b*V(i-1)
! Gram-Schmitt's reorthogonalisation
do j=1,i-1
dotprod= psb_gedot(V(i),V(j),desc_a,info) ! dotprod = (V(i) dot V(j))
call psb_geaxpby(-dotprod,V(j),done,V(i),desc_a, info) !V(i)=V(i)-V(j)*dotprod
H(j,i-1)=cmplx(dotprod,0.0)
end do
norm = psb_norm2(V(i),desc_a,info)
! write(psb_out_unit,'("norma finale :"i20,F20.3)')i,norm
! write(psb_out_unit,'("")')
if (i .ne. dim_H+1) then
H(i,i-1)=cmplx(norm,0.0)
endif
norm=1/norm
call psb_geaxpby(dzero,V(i),norm,V(i),desc_a, info)
enddo
write(psb_out_unit,'("")')
if (iam==psb_root_) then
allocate(eig(dim_H),work(dim_h),Z(dim_H,dim_H),stat = info)
call ZHSEQR('E','N',dim_H,1,dim_H,H,dim_H,eig,Z,dim_H,work,dim_H,info)
!sort H's eigenvalues
allocate(indexes(1:dim_H))
call psb_qsort(eig,indexes,psb_alsort_up_,psb_sort_ovw_idx_)
end if
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
call psb_amx(ictxt,t2)
nr = desc_a%get_global_rows()
annz = a%get_nzeros()
amatsize = psb_sizeof(a)
descsize = psb_sizeof(desc_a)
call psb_sum(ictxt,annz)
call psb_sum(ictxt,amatsize)
call psb_sum(ictxt,descsize)
if (iam==psb_root_) then
flops = 2.d0*times*annz
tflops=flops
! write(psb_out_unit,'("Matrix: ",a)') mtrx_file
! write(psb_out_unit,'("Test on : ",i20," processors")') np
! write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
! write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
open(15, FILE="resultats.dat", position = 'append',ACTION="WRITE")
write (15, '(a,F20.6,F20.6,F20.4)')mtrx_file,real(eig(dim_H)),real(eig(dim_H-1)),t2
close(15)
DEALLOCATE (work,eig,Z)
!write(*,'("valeurs propres de H : ")')
!do i=dim_H/3,dim_H
! write(psb_out_unit,'(g20.4,g20.4)')real(eig(i)),aimag(eig(i))
!enddo
tflops = tflops / (tt2)
!
! This computation is valid for CSR
!
nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_int)+ &
& annz*(psb_sizeof_sp + psb_sizeof_int)
bdwdth = times*nbytes/(t2*1.d6)
bdwdth = times*nbytes/(tt2*1.d6)
end if
call psb_gefree(b_col, desc_a,info)
call psb_gefree(x_col, desc_a,info)
call psb_spfree(a, desc_a,info)
call psb_cdfree(desc_a,info)
do i=1,dim_H
call psb_gefree(V(i), desc_a,info)
enddo
DEALLOCATE (H)
DEALLOCATE (V)
9999 continue
if(info /= 0) then
call psb_error(ictxt)
end if
call psb_exit(ictxt)
stop
contains
subroutine adj_read (a,filename,iunit,desc_a,info)
type(psb_dspmat_type), intent (inout) :: a
character(len=40) :: filename
integer (psb_ipk_) :: iunit
type(psb_desc_type):: desc_a
integer (psb_ipk_) :: info
integer(psb_ipk_) :: i,nnzero,nrows
integer (psb_ipk_) :: iError
type(psb_d_coo_sparse_mat) :: acoo
open(iunit, FILE=filename, STATUS="OLD", ACTION="READ")
read(iunit, *) nrows , nnzero
call acoo%allocate(nrows,nrows,nnzero)
do i = 1,nnzero
read(iunit, *) acoo%ia(i),acoo%ja(i)
acoo%ia(i)=acoo%ia(i)+1
acoo%ja(i)=acoo%ja(i)+1
acoo%val(i)=1.0
end do
close(UNIT=iunit)
!call psb_spall(a,desc_a,info,nnzero)
!call psb_spins(nnzero, ia, ja, val, a, desc_a, info)
call acoo%set_nzeros(nnzero)
call acoo%fix(info)
call a%mv_from(acoo)
call a%cscnv(info,type='csr')
end subroutine adj_read
subroutine lapl(a,b)
type(psb_dspmat_type),intent(in)::a
type(psb_dspmat_type),intent(out)::b
type(psb_d_coo_sparse_mat) :: acoo
integer(psb_ipk_) :: nz,n,info,i
real(psb_dpk_), allocatable :: K(:)
call a%cp_to(acoo)
nz=acoo%get_nzeros()
n=a%get_nrows()
allocate(K(n))
do i=1,n
K(i)=0
enddo
do i=1,nz
K(acoo%ia(i))=K(acoo%ia(i))+acoo%val(i)
acoo%val(i)=-acoo%val(i)
enddo
call acoo%reallocate(nz+n)
do i=1,n
acoo%val(nz+i)=K(i)
acoo%ia(nz+i)= i
acoo%ja(nz+i)= i
enddo
call acoo%set_nzeros(nz+n)
call acoo%fix(info)
call b%mv_from(acoo)
!if(iam==psb_root_) then
! write(psb_out_unit,'("nz, n",i10,i10," 1")') nz,n
! write(psb_out_unit,'("b%get_nzeros ?",i10, " 1")') b%get_nzeros()
!else
! write(psb_out_unit,'("nz,n",i10,i10," 2")') nz,n
! write(psb_out_unit,'("b%get_nzeros ?",i10, " 2")')b%get_nzeros()
!end if
call b%cscnv(info,'CSR')
deallocate (K)
end subroutine lapl
end program arnoldi_file

284
test/kernel/laplacian.f90
Unescape Escape View File

@ -0,0 +1,284 @@
program laplacian
use psb_base_mod
use psb_util_mod
implicit none
! input parameters
character(len=40) :: mtrx_file, rhs_file
! sparse matrices
type(psb_dspmat_type) :: a, b, aux_a
! dense matrices
real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:)
real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:)
real(psb_dpk_), pointer :: b_col_glob(:)
type(psb_d_vect_type) :: b_col, x_col, r_col
! communications data structure
type(psb_desc_type):: desc_a
integer(psb_ipk_) :: ictxt, iam, np
! solver paramters
integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
& methd, istopc, irst, nr
integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes
real(psb_dpk_) :: err, eps,cond
character(len=5) :: afmt
character(len=20) :: name
character(len=2) :: filefmt
integer(psb_ipk_), parameter :: iunit=12
integer(psb_ipk_) :: times=0
integer(psb_ipk_) :: iparm(20)
! other variables
integer(psb_ipk_) :: i,info,j,m_problem
integer(psb_ipk_) :: internal, m,ii,nnzero
real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
&scale,resmx,resmxp, flops, bdwdth
real(psb_dpk_) :: tt1, tt2, tflops
real(psb_dpk_) :: lambda, lambda2
real (psb_dpk_) :: norm, precisione
integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne
integer(psb_ipk_), allocatable :: ivg(:), ipv(:)
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ictxt)
stop
endif
name='laplacian'
if(psb_get_errstatus() /= 0) goto 9999
info=psb_success_
call psb_set_errverbosity(2)
!
! Hello world
!
if (iam == psb_root_) then
!write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
!write(*,*) 'This is the ',trim(name),' sample program'
read(psb_inp_unit,*) mtrx_file
read(psb_inp_unit,*) filefmt
read(psb_inp_unit,*) ipart
read(psb_inp_unit,*) precisione
!write (psb_out_unit, '("The precision of the power method is ",F30.20)') precisione
end if
call psb_bcast(ictxt,mtrx_file)
call psb_bcast(ictxt,filefmt)
call psb_bcast(ictxt,ipart)
call psb_bcast(ictxt,precisione)
rhs_file = 'NONE'
afmt = 'CSR'
call psb_barrier(ictxt)
t1 = psb_wtime()
! read the input matrix to be processed and (possibly) the rhs
nrhs = 1
if (iam==psb_root_) then
select case(psb_toupper(filefmt))
case('MM')
! For Matrix Market we have an input file for the matrix
! and an (optional) second file for the RHS.
call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
if (info == psb_success_) then
if (rhs_file /= 'NONE') then
call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
end if
end if
case ('HB')
! For Harwell-Boeing we have a single file which may or may not
! contain an RHS.
call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
case ('AD')
call adj_read(aux_a,mtrx_file,iunit,desc_a,info)
case default
info = -1
write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
end select
if (info /= psb_success_) then
write(psb_err_unit,*) 'Error while reading input matrix '
call psb_abort(ictxt)
end if
m_problem = aux_a%get_nrows()
call psb_bcast(ictxt,m_problem)
! At this point aux_b may still be unallocated
if (psb_size(aux_b,dim=1)==m_problem) then
! if any rhs were present, broadcast the first one
write(psb_err_unit,'("Ok, got an rhs ")')
b_col_glob =>aux_b(:,1)
else
!write(psb_out_unit,'("Generating an rhs...")')
!write(psb_out_unit,'(" ")')
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob => aux_b(:,1)
do i=1, m_problem
b_col_glob(i) = 1.d0
enddo
endif
call psb_bcast(ictxt,b_col_glob(1:m_problem))
else
call psb_bcast(ictxt,m_problem)
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob =>aux_b(:,1)
call psb_bcast(ictxt,b_col_glob(1:m_problem))
end if
! switch over different partition types
if (ipart == 0) then
call psb_barrier(ictxt)
!if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
allocate(ivg(m_problem),ipv(np))
do i=1,m_problem
call part_block(i,m_problem,np,ipv,nv)
ivg(i) = ipv(1)
enddo
call psb_matdist(aux_a, a,ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else if (ipart == 2) then
if (iam==psb_root_) then
!write(psb_out_unit,'("Partition type: graph")')
!write(psb_out_unit,'(" ")')
! write(psb_err_unit,'("Build type: graph")')
call build_mtpart(aux_a,np)
endif
call psb_barrier(ictxt)
call distr_mtpart(psb_root_,ictxt)
call getv_mtpart(ivg)
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else
!if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block)
end if
call lapl(a,b)
! call psb_gefree(b_col, desc_a,info)
!call psb_gefree(x_col, desc_a,info)
call psb_spfree(a, desc_a,info)
call psb_spfree(b, desc_a,info)
call psb_cdfree(desc_a,info)
9999 continue
if(info /= 0) then
call psb_error(ictxt)
end if
call psb_exit(ictxt)
stop
contains
subroutine adj_read (a,filename,iunit,desc_a,info)
type(psb_dspmat_type), intent (inout) :: a
character(len=40) :: filename
integer (psb_ipk_) :: iunit
type(psb_desc_type):: desc_a
integer (psb_ipk_) :: info
integer(psb_ipk_) :: i,nnzero,nrows
integer (psb_ipk_) :: iError
type(psb_d_coo_sparse_mat) :: acoo
open(iunit, FILE=filename, STATUS="OLD", ACTION="READ")
read(iunit, *) nrows , nnzero
call acoo%allocate(nrows,nrows,nnzero)
do i = 1,nnzero
read(iunit, *) acoo%ia(i),acoo%ja(i)
acoo%ia(i)=acoo%ia(i)+1
acoo%ja(i)=acoo%ja(i)+1
acoo%val(i)=1.0
end do
close(UNIT=iunit)
call acoo%set_nzeros(nnzero)
call acoo%fix(info)
call a%mv_from(acoo)
call a%cscnv(info,type='csr')
end subroutine adj_read
subroutine lapl(a,b)
type(psb_dspmat_type),intent(in)::a
type(psb_dspmat_type),intent(out)::b
type(psb_d_coo_sparse_mat) :: acoo
integer(psb_ipk_) :: nz,n,info,i
real(psb_dpk_), allocatable :: K(:)
call a%cp_to(acoo)
nz=acoo%get_nzeros()
n=a%get_nrows()
allocate(K(n))
do i=1,n
K(i)=0
enddo
do i=1,nz
K(acoo%ia(i))=K(acoo%ia(i))+acoo%val(i)
acoo%val(i)=-acoo%val(i)
enddo
call acoo%reallocate(nz+n)
call acoo%set_dupl(psb_dupl_add_)
do i=1,n
acoo%val(nz+i)=K(i)
acoo%ia(nz+i)= i
acoo%ja(nz+i)= i
enddo
call acoo%set_nzeros(nz+n)
call acoo%fix(info)
do i=1,nz
if(acoo%ja(i)==acoo%ia(i)) then
write(psb_out_unit,'(i10,i10,g20.4)')acoo%ia(i),acoo%ja(i), acoo%val(i)
end if
enddo
!if(iam==psb_root_) then
! write(psb_out_unit,'("nz, n",i10,i10," 1")') nz,n
! write(psb_out_unit,'("b%get_nzeros ?",i10, " 1")')b%get_nzeros()
! else
! write(psb_out_unit,'("nz,n",i10,i10," 2")') nz,n
! write(psb_out_unit,'("b%get_nzeros ?",i10, " 2")')b%get_nzeros()
!end if
call b%mv_from(acoo)
! if(iam==psb_root_) then
! write(psb_out_unit,'("nz, n",i10,i10," 1")') nz,n
! write(psb_out_unit,'("b%get_nzeros ?",i10, " 1")')b%get_nzeros()
! else
! write(psb_out_unit,'("nz,n",i10,i10," 2")') nz,n
! write(psb_out_unit,'("b%get_nzeros ?",i10, " 2")')b%get_nzeros()
! end if
call b%cscnv(info,'CSR')
deallocate (K)
end subroutine lapl
end program laplacian

300
test/kernel/power_file.f90
Unescape Escape View File

@ -0,0 +1,300 @@
!!$
!!$ Parallel Sparse BLAS version 3.1
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
program power_file
use psb_base_mod
use psb_util_mod
implicit none
! input parameters
character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
! sparse matrices
type(psb_zspmat_type) :: a, aux_a
! dense matrices
complex(psb_dpk_), allocatable, target :: aux_b(:,:), d(:)
complex(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:)
complex(psb_dpk_), pointer :: b_col_glob(:)
type(psb_z_vect_type) :: b_col, x_col, r_col
! communications data structure
type(psb_desc_type):: desc_a
integer(psb_ipk_) :: ictxt, iam, np
! solver paramters
integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
& methd, istopc, irst, nr
integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes
real(psb_dpk_) :: err, eps,cond
character(len=5) :: afmt
character(len=20) :: name
character(len=2) :: filefmt
integer(psb_ipk_), parameter :: iunit=12
integer(psb_ipk_) :: times=0
integer(psb_ipk_) :: iparm(20)
! other variables
integer(psb_ipk_) :: i,info,j,m_problem
integer(psb_ipk_) :: internal, m,ii,nnzero
real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
&scale,resmx,resmxp, flops, bdwdth
real(psb_dpk_) :: tt1, tt2, tflops
complex(psb_dpk_) :: lambda, lambda2
real (psb_dpk_) :: norm, precisione
integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne
integer(psb_ipk_), allocatable :: ivg(:), ipv(:)
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ictxt)
stop
endif
name='power_file'
if(psb_get_errstatus() /= 0) goto 9999
info=psb_success_
call psb_set_errverbosity(2)
!
! Hello world
!
if (iam == psb_root_) then
write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
write(*,*) 'This is the ',trim(name),' sample program'
read(psb_inp_unit,*) mtrx_file
read(psb_inp_unit,*) filefmt
read(psb_inp_unit,*) ipart
read(psb_inp_unit,*) precisione
write (psb_out_unit, '("The precision of the power method is ",F30.20)') precisione
end if
call psb_bcast(ictxt,mtrx_file)
call psb_bcast(ictxt,filefmt)
call psb_bcast(ictxt,ipart)
call psb_bcast(ictxt,precisione)
rhs_file = 'NONE'
afmt = 'CSR'
call psb_barrier(ictxt)
t1 = psb_wtime()
! read the input matrix to be processed and (possibly) the rhs
nrhs = 1
if (iam==psb_root_) then
select case(psb_toupper(filefmt))
case('MM')
! For Matrix Market we have an input file for the matrix
! and an (optional) second file for the RHS.
call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
if (info == psb_success_) then
if (rhs_file /= 'NONE') then
call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
end if
end if
case ('HB')
! For Harwell-Boeing we have a single file which may or may not
! contain an RHS.
call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
case default
info = -1
write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
end select
if (info /= psb_success_) then
write(psb_err_unit,*) 'Error while reading input matrix '
call psb_abort(ictxt)
end if
m_problem = aux_a%get_nrows()
call psb_bcast(ictxt,m_problem)
! At this point aux_b may still be unallocated
if (psb_size(aux_b,dim=1)==m_problem) then
! if any rhs were present, broadcast the first one
write(psb_err_unit,'("Ok, got an rhs ")')
b_col_glob =>aux_b(:,1)
else
write(psb_out_unit,'("Generating an rhs...")')
write(psb_out_unit,'(" ")')
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob => aux_b(:,1)
do i=1, m_problem
b_col_glob(i) = 1.d0
enddo
endif
call psb_bcast(ictxt,b_col_glob(1:m_problem))
else
call psb_bcast(ictxt,m_problem)
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob =>aux_b(:,1)
call psb_bcast(ictxt,b_col_glob(1:m_problem))
end if
! switch over different partition types
if (ipart == 0) then
call psb_barrier(ictxt)
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
allocate(ivg(m_problem),ipv(np))
do i=1,m_problem
call part_block(i,m_problem,np,ipv,nv)
ivg(i) = ipv(1)
enddo
call psb_matdist(aux_a, a,ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else if (ipart == 2) then
if (iam==psb_root_) then
write(psb_out_unit,'("Partition type: graph")')
write(psb_out_unit,'(" ")')
! write(psb_err_unit,'("Build type: graph")')
call build_mtpart(aux_a,np)
endif
call psb_barrier(ictxt)
call distr_mtpart(psb_root_,ictxt)
call getv_mtpart(ivg)
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block)
end if
call psb_geall(x_col,desc_a,info)
call x_col%set(zone)
call psb_geasb(x_col,desc_a,info)
t2 = psb_wtime() - t1
call psb_amx(ictxt, t2)
if (iam==psb_root_) then
write(psb_out_unit,'(" ")')
write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
write(psb_out_unit,'(" ")')
end if
call psb_barrier(ictxt)
t1 = psb_wtime()
call psb_spmm(zone,a,x_col,zzero,b_col,desc_a,info,'n') !we do b_col=a*x_col
norm = psb_norm2(b_col,desc_a,info)
norm = 1/norm
!normalisation of b_col in x_col
call psb_geaxpby(zzero,x_col,norm*zone,b_col,desc_a, info)
x_col=b_col
lambda=0
lambda2=1000
do while (abs(lambda-lambda2)>precisione)
!write (psb_out_unit, '("times",i20)')times
times=times+1
lambda2=lambda
call psb_spmm(zone,a,x_col,zzero,b_col,desc_a,info,'n')
norm = psb_norm2(b_col,desc_a,info)
lambda = psb_gedot(b_col,x_col,desc_a,info) ! lambda = (A*x_col dot x_col)
norm = 1/norm
call psb_geaxpby(zzero,x_col,norm*zone,b_col,desc_a, info)
x_col=b_col
!write (psb_out_unit, '("abs(lambda-lambda2)",F30.20)') abs(lambda-lambda2)
!write (psb_out_unit, '("lambda",F30.20,F30.20)') lambda
! write (psb_out_unit, '("precisione",F30.20,i20)') precisione, times
end do
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
call psb_amx(ictxt,t2)
nr = desc_a%get_global_rows()
annz = a%get_nzeros()
amatsize = psb_sizeof(a)
descsize = psb_sizeof(desc_a)
call psb_sum(ictxt,annz)
call psb_sum(ictxt,amatsize)
call psb_sum(ictxt,descsize)
if (iam==psb_root_) then
flops = 2.d0*times*annz
tflops=flops
write(psb_out_unit,'("Matrix: ",a)') mtrx_file
write(psb_out_unit,'("Test on : ",i20," processors")') np
write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize
write(psb_out_unit,'("Number of flops (",i0," iters) : ",F20.0," ")') times,flops
flops = flops / (t2)
tflops = tflops / (tt2)
write (psb_out_unit, '("The biggest eigenvalue of A is : ",F20.3,F20.3)')lambda
!
! This computation is valid for CSR
!
nbytes = nr*(2*psb_sizeof_sp + psb_sizeof_int)+ &
& annz*(psb_sizeof_sp + psb_sizeof_int)
bdwdth = times*nbytes/(t2*1.d6)
bdwdth = times*nbytes/(tt2*1.d6)
end if
call psb_gefree(b_col, desc_a,info)
call psb_gefree(x_col, desc_a,info)
call psb_spfree(a, desc_a,info)
call psb_cdfree(desc_a,info)
9999 continue
if(info /= 0) then
call psb_error(ictxt)
end if
call psb_exit(ictxt)
stop
end program power_file

343
test/kernel/power_file_real.f90
Unescape Escape View File

@ -0,0 +1,343 @@
!!$
!!$ Parallel Sparse BLAS version 3.1
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
program power_file_real
use psb_base_mod
use psb_util_mod
implicit none
! input parameters
character(len=100) :: kmethd, ptype, mtrx_file, rhs_file
! sparse matrices
type(psb_dspmat_type) :: a, aux_a
! dense matrices
real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:)
real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:)
real(psb_dpk_), pointer :: b_col_glob(:)
type(psb_d_vect_type) :: b_col, x_col, r_col
! communications data structure
type(psb_desc_type):: desc_a
integer(psb_ipk_) :: ictxt, iam, np
! solver paramters
integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
& methd, istopc, irst, nr
integer(psb_long_int_k_) :: amatsize, descsize, annz, nbytes
real(psb_dpk_) :: err, eps,cond
character(len=5) :: afmt
character(len=20) :: name
character(len=2) :: filefmt
integer(psb_ipk_), parameter :: iunit=12
integer(psb_ipk_) :: times=0
integer(psb_ipk_) :: iparm(20)
! other variables
integer(psb_ipk_) :: i,info,j,m_problem
integer(psb_ipk_) :: internal, m,ii,nnzero
real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
&scale,resmx,resmxp, flops, bdwdth
real(psb_dpk_) :: tt1, tt2, tflops
real(psb_dpk_) :: lambda, lambda2
real (psb_dpk_) :: norm, precisione
integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne
integer(psb_ipk_), allocatable :: ivg(:), ipv(:)
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ictxt)
stop
endif
name='power_file'
if(psb_get_errstatus() /= 0) goto 9999
info=psb_success_
call psb_set_errverbosity(2)
!
! Hello world
!
if (iam == psb_root_) then
!write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
!write(*,*) 'This is the ',trim(name),' sample program'
read(psb_inp_unit,*) mtrx_file
read(psb_inp_unit,*) filefmt
read(psb_inp_unit,*) ipart
read(psb_inp_unit,*) precisione
!write (psb_out_unit, '("The precision of the power method is ",F30.20)') precisione
end if
call psb_bcast(ictxt,mtrx_file)
call psb_bcast(ictxt,filefmt)
call psb_bcast(ictxt,ipart)
call psb_bcast(ictxt,precisione)
rhs_file = 'NONE'
afmt = 'CSR'
call psb_barrier(ictxt)
t1 = psb_wtime()
! read the input matrix to be processed and (possibly) the rhs
nrhs = 1
if (iam==psb_root_) then
select case(psb_toupper(filefmt))
case('MM')
! For Matrix Market we have an input file for the matrix
! and an (optional) second file for the RHS.
call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
if (info == psb_success_) then
if (rhs_file /= 'NONE') then
call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
end if
end if
case ('HB')
! For Harwell-Boeing we have a single file which may or may not
! contain an RHS.
call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
case ('AD')
call adj_read(aux_a,mtrx_file,iunit,desc_a,info)
case default
info = -1
write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
end select
if (info /= psb_success_) then
write(psb_err_unit,*) 'Error while reading input matrix '
call psb_abort(ictxt)
end if
m_problem = aux_a%get_nrows()
call psb_bcast(ictxt,m_problem)
! At this point aux_b may still be unallocated
if (psb_size(aux_b,dim=1)==m_problem) then
! if any rhs were present, broadcast the first one
write(psb_err_unit,'("Ok, got an rhs ")')
b_col_glob =>aux_b(:,1)
else
!write(psb_out_unit,'("Generating an rhs...")')
!write(psb_out_unit,'(" ")')
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob => aux_b(:,1)
do i=1, m_problem
b_col_glob(i) = 1.d0
enddo
endif
call psb_bcast(ictxt,b_col_glob(1:m_problem))
else
call psb_bcast(ictxt,m_problem)
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob =>aux_b(:,1)
call psb_bcast(ictxt,b_col_glob(1:m_problem))
end if
! switch over different partition types
if (ipart == 0) then
call psb_barrier(ictxt)
!if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
allocate(ivg(m_problem),ipv(np))
do i=1,m_problem
call part_block(i,m_problem,np,ipv,nv)
ivg(i) = ipv(1)
enddo
call psb_matdist(aux_a, a,ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else if (ipart == 2) then
if (iam==psb_root_) then
!write(psb_out_unit,'("Partition type: graph")')
!write(psb_out_unit,'(" ")')
! write(psb_err_unit,'("Build type: graph")')
call build_mtpart(aux_a,np)
endif
call psb_barrier(ictxt)
call distr_mtpart(psb_root_,ictxt)
call getv_mtpart(ivg)
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else
!if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block)
end if
call psb_geall(x_col,desc_a,info)
call x_col%set(done)
call psb_geasb(x_col,desc_a,info)
t2 = psb_wtime() - t1
call psb_amx(ictxt, t2)
if (iam==psb_root_) then
!write(psb_out_unit,'(" ")')
!write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
!write(psb_out_unit,'(" ")')
end if
call psb_barrier(ictxt)
call psb_spmm(done,a,x_col,dzero,b_col,desc_a,info,'n') !we do b_col=a*x_col
norm = psb_norm2(b_col,desc_a,info)
norm = 1/norm
!normalisation of b_col in x_col
call psb_geaxpby(dzero,x_col,norm,b_col,desc_a, info)
x_col=b_col
lambda=0
lambda2=1000
do while (abs(lambda-lambda2)>precisione .AND. times<300)
!if(iam==psb_root_) write (psb_out_unit,'(F20.6,i20)')abs(lambda-lambda2), times
times=times+1
lambda2=lambda
call psb_spmm(done,a,x_col,dzero,b_col,desc_a,info,'n')
norm = psb_norm2(b_col,desc_a,info)
lambda = psb_gedot(b_col,x_col,desc_a,info) ! lambda = (A*x_col dot x_col)
norm = 1/norm
call psb_geaxpby(dzero,x_col,norm,b_col,desc_a, info)
x_col=b_col
end do
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
call psb_amx(ictxt,t2)
nr = desc_a%get_global_rows()
annz = a%get_nzeros()
amatsize = psb_sizeof(a)
descsize = psb_sizeof(desc_a)
call psb_sum(ictxt,annz)
call psb_sum(ictxt,amatsize)
call psb_sum(ictxt,descsize)
if (iam==psb_root_) then
flops = 2.d0*times*annz
tflops=flops
! write(psb_out_unit,'("Matrix: ",a)') mtrx_file
! write(psb_out_unit,'("Test on : ",i20," processors")') np
! write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
! write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
!write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize
! write(psb_out_unit,'("Number of flops (",i0," iters) : ",F20.0," ")') times,flops
! write(psb_out_unit,'("Time for ",i0," iteration (s) : ",F20.3)')times, t2
! write(psb_out_unit,'("Time per iteration (ms) : ",F20.3)') t2*1.d3/(1.d0*times)
open(15, FILE="resultats.dat", position = 'append',ACTION="WRITE")
write (15, '(a,F20.6,F20.4)')mtrx_file,lambda,t2
! write (psb_out_unit, '("entropy vect propre :",g20.4)')entropy(x_col)
close(15)
end if
call psb_gefree(b_col, desc_a,info)
call psb_gefree(x_col, desc_a,info)
call psb_spfree(a, desc_a,info)
call psb_cdfree(desc_a,info)
9999 continue
if(info /= 0) then
call psb_error(ictxt)
end if
call psb_exit(ictxt)
stop
contains
subroutine adj_read (a,filename,iunit,desc_a,info)
type(psb_dspmat_type), intent (inout) :: a
character(len=100) :: filename
integer (psb_ipk_) :: iunit
type(psb_desc_type):: desc_a
integer (psb_ipk_) :: info
integer(psb_ipk_) :: i,nnzero,nrows
integer (psb_ipk_) :: iError
!integer(psb_ipk_), allocatable :: ia(:),ja(:)
!real (psb_dpk_), allocatable :: val(:)
type(psb_d_coo_sparse_mat) :: acoo
open(iunit, FILE=filename, STATUS="OLD", ACTION="READ")
read(iunit, *) nrows , nnzero
call acoo%allocate(nrows,nrows,nnzero)
do i = 1,nnzero
read(iunit, *) acoo%ia(i),acoo%ja(i)
acoo%ia(i)=acoo%ia(i)+1
acoo%ja(i)=acoo%ja(i)+1
acoo%val(i)=1.0
end do
close(UNIT=iunit)
!call psb_spall(a,desc_a,info,nnzero)
!call psb_spins(nnzero, ia, ja, val, a, desc_a, info)
call acoo%set_nzeros(nnzero)
call acoo%fix(info)
call a%mv_from(acoo)
call a%cscnv(info,type='csr')
end subroutine adj_read
real(psb_dpk_) function entropy(vect)
type(psb_d_vect_type), intent(inout)::vect
real(psb_dpk_), allocatable :: v(:)
integer(psb_ipk_) :: nrows
nrows=vect%get_nrows()
allocate(v(nrows))
v=vect%get_vect()
entropy=0
do i=1,nrows
entropy=entropy+v(i)*v(i)*log(v(i)*v(i))
end do
return
end function entropy
end program power_file_real

471
test/kernel/shift_invert_real.f90
Unescape Escape View File

@ -0,0 +1,471 @@
!!$
!!$ Parallel Sparse BLAS version 3.1
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
program shift_invert
use psb_base_mod
use psb_prec_mod
use psb_krylov_mod
use psb_util_mod
implicit none
! input parameters
character(len=40) :: kmethd, ptype, mtrx_file, rhs_file
! sparse matrices and preconditioner
type(psb_dspmat_type) :: a, aux_a, id, a2, b
type(psb_d_csr_sparse_mat)::acsr, icsr
type(psb_dprec_type) :: prec
! dense matrices
real(psb_dpk_), allocatable, target :: aux_b(:,:), d(:)
complex(psb_dpk_), allocatable, target :: H(:,:),eig(:),work(:),Z(:,:)
real(psb_dpk_), allocatable , save :: x_col_glob(:), r_col_glob(:)
real(psb_dpk_), pointer :: b_col_glob(:)
type(psb_d_vect_type) :: b_col, x_col, r_col
type (psb_d_vect_type), allocatable, target :: V(:)
integer(psb_ipk_), allocatable , save :: indexes(:)
! communications data structure
type(psb_desc_type):: desc_a
integer(psb_ipk_) :: ictxt, iam, np
! solver paramters
integer(psb_ipk_) :: iter, itmax, ierr, itrace, ircode, ipart,&
& methd, istopc, irst, nr, sort, annz
integer(psb_long_int_k_) :: amatsize, descsize, nbytes
real(psb_dpk_) :: err, eps,cond
character(len=5) :: afmt
character(len=20) :: name, ch_err
character(len=2) :: filefmt
integer(psb_ipk_), parameter :: iunit=12
integer(psb_ipk_) :: times=0
integer(psb_ipk_) :: iparm(20)
! other variables
integer(psb_ipk_) :: i,info,j,m_problem
integer(psb_ipk_) :: internal, m,ii,nnzero,ji
real(psb_dpk_) :: t1, t2, r_amax, b_amax,&
&scale,resmx,resmxp, flops, bdwdth
real(psb_dpk_) :: tt1, tt2, tflops, tprec, sigma,dotprod,norm
integer(psb_ipk_) :: nrhs, nrow, n_row, dim, nv, ne,dim_H
integer(psb_ipk_), allocatable :: ivg(:), ipv(:)
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ictxt)
stop
endif
name='shift_invert_real'
if(psb_get_errstatus() /= 0) goto 9999
info=psb_success_
call psb_set_errverbosity(2)
!
! Hello world
!
if (iam == psb_root_) then
write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
write(*,*) 'This is the ',trim(name),' sample program'
read(psb_inp_unit,*) mtrx_file
read(psb_inp_unit,*) filefmt
read(psb_inp_unit,*) ipart
read(psb_inp_unit,*) ptype
read(psb_inp_unit,*) kmethd
read(psb_inp_unit,*) dim_H
read(psb_inp_unit,*) sigma
read(psb_inp_unit,*) eps
end if
call psb_bcast(ictxt,mtrx_file)
call psb_bcast(ictxt,filefmt)
call psb_bcast(ictxt,ipart)
call psb_bcast(ictxt,dim_H)
call psb_bcast(ictxt,kmethd)
call psb_bcast(ictxt,ptype)
call psb_bcast(ictxt,sigma)
call psb_bcast(ictxt,eps)
rhs_file = 'NONE'
afmt = 'CSR'
call psb_barrier(ictxt)
t1 = psb_wtime()
! read the input matrix to be processed and (possibly) the rhs
nrhs = 1
if (iam==psb_root_) then
select case(psb_toupper(filefmt))
case('MM')
! For Matrix Market we have an input file for the matrix
! and an (optional) second file for the RHS.
call mm_mat_read(aux_a,info,iunit=iunit,filename=mtrx_file)
if (info == psb_success_) then
if (rhs_file /= 'NONE') then
call mm_array_read(aux_b,info,iunit=iunit,filename=rhs_file)
end if
end if
case ('HB')
! For Harwell-Boeig we have a single file which may or may not
! contain an RHS.
call hb_read(aux_a,info,iunit=iunit,b=aux_b,filename=mtrx_file)
case ('AD')
call adj_read(aux_a,mtrx_file,iunit,desc_a,info)
case default
info = -1
write(psb_err_unit,*) 'Wrong choice for fileformat ', filefmt
end select
if (info /= psb_success_) then
write(psb_err_unit,*) 'Error while reading input matrix '
call psb_abort(ictxt)
end if
m_problem = aux_a%get_nrows()
annz=aux_a%get_nzeros()
call psb_bcast(ictxt,m_problem)
! At this point aux_b may still be unallocated
if (psb_size(aux_b,dim=1)==m_problem) then
! if any rhs were present, broadcast the first one
write(psb_err_unit,'("Ok, got an rhs ")')
b_col_glob =>aux_b(:,1)
else
! write(psb_out_unit,'("Generating an rhs...")')
! write(psb_out_unit,'(" ")')
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob => aux_b(:,1)
do i=1, m_problem
b_col_glob(i) = 1.d0
enddo
endif
call psb_bcast(ictxt,b_col_glob(1:m_problem))
else
call psb_bcast(ictxt,m_problem)
call psb_realloc(m_problem,1,aux_b,ircode)
if (ircode /= 0) then
call psb_errpush(psb_err_alloc_dealloc_,name)
goto 9999
endif
b_col_glob =>aux_b(:,1)
call psb_bcast(ictxt,b_col_glob(1:m_problem))
end if
! switch over different partition types
if (ipart == 0) then
call psb_barrier(ictxt)
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
allocate(ivg(m_problem),ipv(np))
do i=1,m_problem
call part_block(i,m_problem,np,ipv,nv)
ivg(i) = ipv(1)
enddo
call psb_matdist(aux_a, a,ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else if (ipart == 2) then
if (iam==psb_root_) then
!write(psb_out_unit,'("Partition type: graph")')
!write(psb_out_unit,'(" ")')
! write(psb_err_unit,'("Build type: graph")')
call build_mtpart(aux_a,np)
endif
call psb_barrier(ictxt)
call distr_mtpart(psb_root_,ictxt)
call getv_mtpart(ivg)
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,v=ivg)
else
if (iam==psb_root_) write(psb_out_unit,'("Partition type: block")')
call psb_matdist(aux_a, a, ictxt, &
& desc_a,b_col_glob,b_col,info,fmt=afmt,parts=part_block)
end if
call lapl(a,b)
!
!a2=a-sigma*id
!
!sigma = 100
call csshift(b,a2,sigma,ictxt)
!
! prepare the preconditioner.
!
!if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ", a) ') ptype
call psb_precinit(prec,ptype,info)
call psb_barrier(ictxt)
!t1 = psb_wtime()
call psb_precbld(a2,desc_a,prec,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_precbld'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
tprec = psb_wtime()-t1
if (iam == psb_root_) write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec
allocate(H(dim_H,dim_H))
allocate (V(dim_H+1))
do i=1,dim_H+1
call psb_geall(V(i),desc_a,info)
call psb_geasb(V(i),desc_a,info)
enddo
call V(1)%set(done)
call psb_amx(ictxt, t2)
!if (iam==psb_root_) then
! write(psb_out_unit,'(" ")')
! write(psb_out_unit,'("Time to read and partition matrix : ",es12.5)')t2
! write(psb_out_unit,'(" ")')
!end if
call psb_barrier(ictxt)
t2 = psb_wtime()-t1-tprec
if (iam == psb_root_) write(psb_out_unit,'("Preconditioner time : ",es12.5)')t2
norm = psb_norm2(V(1),desc_a, info)
H(2,1)=cmplx(norm,0.0)
norm = 1/norm
call psb_geaxpby(dzero,V(1),norm,V(1),desc_a,info)
do i=2,dim_H+1
!A*V(i)=V(i-1)
call psb_krylov(kmethd,a2,prec,V(i-1),V(i),eps,desc_a,info,&
& itmax=1000,iter=iter,err=err,itrace=-1,istop=2,irst=002)
!if (iam==psb_root_) write(*,'("iter : "i20)') iter
! Gram-Schmitt's reorthogonalisation
do j=1,i-1
dotprod= psb_gedot(V(i),V(j),desc_a,info) ! dotprod = (V(i) dot V(j))
call psb_geaxpby(-dotprod,V(j),done,V(i),desc_a, info)!V(i)=V(i)-V(j)*dotprod
H(j,i-1)=cmplx(dotprod,0.0)
end do
norm = psb_norm2(V(i),desc_a,info)
if (i .ne. dim_H+1) then
H(i,i-1)=cmplx(norm,0.0)
endif
norm=1/norm
call psb_geaxpby(dzero,V(i),norm,V(i),desc_a, info)
enddo
t2=psb_wtime()-t1-tprec
if (iam==psb_root_) write (psb_out_unit,'("temps de arnoldi : " ,es12.5 )') t2
if(iam==psb_root_) then
allocate(eig(dim_H),work(dim_h),Z(dim_H,dim_H),stat = info)
call ZHSEQR('E','N',dim_H,1,dim_H,H,dim_H,eig,Z,dim_H,work,dim_H,info)
!sort H's eigenvalues
allocate(indexes(1:dim_H))
call psb_qsort(eig,indexes,psb_alsort_up_,psb_sort_ovw_idx_)
do i=1,dim_H
eig (i) = cmplx(sigma,0.0)+1/eig(i)
!write(psb_out_unit, '("eig(i), i", g20.4, i10)')real(eig(i)),i
enddo
end if
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
call psb_amx(ictxt,t2)
nr = desc_a%get_global_rows()
amatsize = psb_sizeof(a)
descsize = psb_sizeof(desc_a)
call psb_sum(ictxt,annz)
call psb_sum(ictxt,amatsize)
call psb_sum(ictxt,descsize)
if (iam==psb_root_) then
flops = 2.d0*times*annz
tflops=flops
!write(psb_out_unit,'("Matrix: ",a)') mtrx_file
!write(psb_out_unit,'("Test on : ",i20," processors")') np
!write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
!write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
!write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize
!write(psb_out_unit,'("Number of flops (",i0," iters) : ",F20.0," ")') times,flops
!write(*,'("eigenvalues near from ", g20.4," : ")') sigma
!do i=dim_H/3,dim_H
! write(psb_out_unit,'(g20.4,g20.4)')real(eig(i)),aimag(eig(i))
!enddo
open(15, FILE="resultats.dat", position = 'append',ACTION="WRITE")
write (15,'(F20.6,F20.6,F20.4)')real(eig(dim_H-1)),real(eig(dim_H)),t2
close(15)
DEALLOCATE (work,eig,Z)
end if
call psb_gefree(b_col, desc_a,info)
call psb_gefree(x_col, desc_a,info)
call psb_spfree(a, desc_a,info)
call psb_spfree(a2, desc_a,info)
call psb_spfree(b, desc_a,info)
call psb_cdfree(desc_a,info)
do i=1,dim_H
call psb_gefree(V(i), desc_a,info)
enddo
DEALLOCATE (H)
DEALLOCATE (V)
9999 continue
if(info /= 0) then
call psb_error(ictxt)
end if
call psb_exit(ictxt)
stop
contains
subroutine csshift(a,b,sigma,ictx)
type(psb_dspmat_type), intent(in) :: a
type(psb_dspmat_type), intent(out) :: b
real(psb_dpk_) :: sigma
integer(psb_ipk_)::ictx
type(psb_d_coo_sparse_mat) :: acoo
integer(psb_ipk_) :: nz,n,info,i
call a%cp_to(acoo)
if (sigma/=0.0) then
nz=acoo%get_nzeros()
n=a%get_nrows()
call acoo%reallocate(nz+n)
call acoo%set_dupl(psb_dupl_add_)
do i=1,n
acoo%val(nz+i)=-sigma
acoo%ia(nz+i)= i
acoo%ja(nz+i)= i
enddo
call acoo%set_nzeros(nz+n)
call acoo%fix(info)
! do i=1,nz
! if(acoo%ja(i)==acoo%ia(i)) then
! write(psb_out_unit,'(i10,i10,g20.4)')acoo%ja(i),i, acoo%val(i)
! end if
! enddo
!write(psb_out_unit,'("autant de nzeros apres fix ?",i10)') acoo%get_nzeros()-nz
end if
call b%mv_from(acoo)
call b%cscnv(info,'CSR')
end subroutine csshift
subroutine adj_read (a,filename,iunit,desc_a,info)
type(psb_dspmat_type), intent (inout) :: a
character(len=40) :: filename
integer (psb_ipk_) :: iunit
type(psb_desc_type):: desc_a
integer (psb_ipk_) :: info
integer(psb_ipk_) :: i,nnzero,nrows
integer (psb_ipk_) :: iError
type(psb_d_coo_sparse_mat) :: acoo
open(iunit, FILE=filename, STATUS="OLD", ACTION="READ")
read(iunit, *) nrows , nnzero
call acoo%allocate(nrows,nrows,nnzero)
do i = 1,nnzero
read(iunit, *) acoo%ia(i),acoo%ja(i)
acoo%ia(i)=acoo%ia(i)+1
acoo%ja(i)=acoo%ja(i)+1
acoo%val(i)=1.0
end do
close(UNIT=iunit)
!call psb_spall(a,desc_a,info,nnzero)
!call psb_spins(nnzero, ia, ja, val, a, desc_a, info)
call acoo%set_nzeros(nnzero)
call acoo%fix(info)
call a%mv_from(acoo)
call a%cscnv(info,type='csr')
end subroutine adj_read
subroutine lapl(a,b)
type(psb_dspmat_type),intent(in)::a
type(psb_dspmat_type),intent(out)::b
type(psb_d_coo_sparse_mat) :: acoo
integer(psb_ipk_) :: nz,n,info,i
real(psb_dpk_), allocatable :: K(:)
call a%cp_to(acoo)
nz=acoo%get_nzeros()
n=a%get_nrows()
allocate(K(n))
do i=1,n
K(i)=0
enddo
do i=1,nz
K(acoo%ia(i))=K(acoo%ia(i))+acoo%val(i)
acoo%val(i)=-acoo%val(i)
enddo
call acoo%reallocate(nz+n)
call acoo%set_dupl(psb_dupl_add_)
do i=1,n
acoo%val(nz+i)=K(i)
acoo%ia(nz+i)= i
acoo%ja(nz+i)= i
enddo
call acoo%set_nzeros(nz+n)
call acoo%fix(info)
do i=1,nz
! if(acoo%ja(i)==acoo%ia(i)) then
! write(psb_out_unit,'(i10,i10,g20.4)')acoo%ia(i),acoo%ja(i),acoo%val(i)
! end if
enddo
call b%mv_from(acoo)
call b%cscnv(info,'CSR')
deallocate (K)
end subroutine lapl
end program shift_invert

20
test/kernel/test_chseqr.f90
Unescape Escape View File

@ -0,0 +1,20 @@
program test_chseqr
complex , dimension (1:3,1:3) :: H,Z,Vr
complex , dimension (1:6,1:6) :: Rwork
complex, dimension (1:3) :: eing,work
integer :: info, N
N=3
H = reshape((/ (1.0,0),(0,0), (2,0), (2,0) , (3,0) , (-4,0) , (0,0) , (0,0) , (2,0) /), shape(H))
do i=1,N
do j=1,N
write(*,'("H : "i5,i5,2x,g20.4,g20.4)')i,j,real(H(i,j)),aimag(H(i,j))
end do
end do
call CHSEQR('E','N',N,1,N,H,N,eing,Z,N,work,N,info)
!call CGEEV('N','N',N,H,N,eing,Z,N,Vr,N,work,3*N,Rwork,info)
do i=1,N
write(*,'("valeur propre de H : "i5, i5, g20.4,g20.4)')info, i,real(eing(i)),aimag(eing(i))
enddo
end program test_chseqr
Write Preview
Loading…
Cancel
Save