1! Copyright (C) 2005 The Scalable Software Infrastructure Project. All rights reserved. 2! 3! Redistribution and use in source and binary forms, with or without 4! modification, are permitted provided that the following conditions are met: 5! 1. Redistributions of source code must retain the above copyright 6! notice, this list of conditions and the following disclaimer. 7! 2. Redistributions in binary form must reproduce the above copyright 8! notice, this list of conditions and the following disclaimer in the 9! documentation and/or other materials provided with the distribution. 10! 3. Neither the name of the project nor the names of its contributors 11! may be used to endorse or promote products derived from this software 12! without specific prior written permission. 13! 14! THIS SOFTWARE IS PROVIDED BY THE SCALABLE SOFTWARE INFRASTRUCTURE PROJECT 15! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 16! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 17! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE SCALABLE SOFTWARE INFRASTRUCTURE 18! PROJECT BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 19! OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 20! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 21! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 22! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 23! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 24! POSSIBILITY OF SUCH DAMAGE. 25 26 implicit none 27 28#include "lisf.h" 29 30 LIS_MATRIX :: A,A0 31 LIS_VECTOR :: b,x,u 32 LIS_SOLVER :: solver 33 LIS_INTEGER :: ia,ierr 34 integer*4 :: nprocs,my_rank 35 LIS_INTEGER :: matrix_type,comm 36 LIS_INTEGER :: omp_get_num_procs,omp_get_max_threads 37 LIS_INTEGER :: m,n,nn,nnz,innz 38 LIS_INTEGER :: i,j,ii,jj,ctr 39 LIS_INTEGER :: is,ie 40 LIS_INTEGER,allocatable :: ptr(:),index(:) 41 LIS_SCALAR,allocatable :: value(:) 42 LIS_SCALAR :: one=1.0d0 43 LIS_INTEGER :: nsol,iter,iter_double,iter_quad 44 real*8 :: time,itime,ptime,p_c_time,p_i_time 45 LIS_REAL :: resid 46 character*256 :: solname,resname 47 character*20 :: m_char,n_char,mt_char,solvername 48 integer*4 :: iargc 49 50 call lis_initialize(ierr) 51 52 comm = LIS_COMM_WORLD 53 54#ifdef USE_MPI 55 call MPI_Comm_size(comm,nprocs,ierr) 56 call MPI_Comm_rank(comm,my_rank,ierr) 57#else 58 nprocs = 1 59 my_rank = 0 60#endif 61 62 matrix_type = LIS_MATRIX_CSR 63 64 ia = iargc() 65 if( ia.lt.5 ) then 66 if( my_rank.eq.0 ) then 67 write(*,'(a)') 'Usage: test2f m n matrix_type solution_filename residual_filename [options]' 68 call lis_finalize(ierr) 69 endif 70 stop 71 endif 72 73 if ( my_rank .eq. 0 ) then 74 write(*,'(a)') '' 75 write(*,'(a,i0)') 'number of processes = ',nprocs 76#ifdef _OPENMP 77 write(*,'(a,i0)') 'max number of threads = ',omp_get_num_procs() 78 write(*,'(a,i0)') 'number of threads = ', omp_get_max_threads() 79#endif 80 endif 81 82! read matrix and vectors from file 83 call getarg(1,m_char) 84 read(m_char,*) m 85 call getarg(2,n_char) 86 read(n_char,*) n 87 call getarg(3,mt_char) 88 read(mt_char,*) matrix_type 89 90 nn = m * n 91 call lis_matrix_create(comm,A,ierr) 92 call CHKERR(ierr) 93 call lis_matrix_set_size(A,0,nn,ierr) 94 call CHKERR(ierr) 95 96 allocate(ptr(0:nn)) 97 allocate(index(0:5*nn-1)) 98 allocate(value(0:5*nn-1)) 99 100 call lis_matrix_get_range(A,is,ie,ierr) 101 ctr = 0 102#ifdef COMPLEX 103 do ii=is-1,ie-2 104 i = ii/m 105 j = ii - i*m 106 if ( i .GT. 0 ) then 107 jj = ii - m 108 index(ctr) = jj 109 value(ctr) = (-1.0d0,0.0d0) 110 ctr = ctr + 1 111 end if 112 if ( i .LT. n-1 ) then 113 jj = ii + m 114 index(ctr) = jj 115 value(ctr) = (-1.0d0,0.0d0) 116 ctr = ctr + 1 117 end if 118 if ( j .GT. 0 ) then 119 jj = ii - 1 120 index(ctr) = jj 121 value(ctr) = (-1.0d0,0.0d0) 122 ctr = ctr + 1 123 end if 124 if ( j .LT. m-1 ) then 125 jj = ii + 1 126 index(ctr) = jj 127 value(ctr) = (-1.0d0,0.0d0) 128 ctr = ctr + 1 129 end if 130 index(ctr) = ii 131 value(ctr) = (4.0d0,0.0d0) 132 ctr = ctr + 1 133 ptr(ii-(is-1)+1) = ctr 134 end do 135#else 136 do ii=is-1,ie-2 137 i = ii/m 138 j = ii - i*m 139 if ( i .GT. 0 ) then 140 jj = ii - m 141 index(ctr) = jj 142 value(ctr) = -1.0d0 143 ctr = ctr + 1 144 end if 145 if ( i .LT. n-1 ) then 146 jj = ii + m 147 index(ctr) = jj 148 value(ctr) = -1.0d0 149 ctr = ctr + 1 150 end if 151 if ( j .GT. 0 ) then 152 jj = ii - 1 153 index(ctr) = jj 154 value(ctr) = -1.0d0 155 ctr = ctr + 1 156 end if 157 if ( j .LT. m-1 ) then 158 jj = ii + 1 159 index(ctr) = jj 160 value(ctr) = -1.0d0 161 ctr = ctr + 1 162 end if 163 index(ctr) = ii 164 value(ctr) = 4.0d0 165 ctr = ctr + 1 166 ptr(ii-(is-1)+1) = ctr 167 end do 168#endif 169 ptr(0) = 0 170 call lis_matrix_set_csr(ptr(ie-is),ptr,index,value,A,ierr) 171 call CHKERR(ierr) 172 call lis_matrix_assemble(A,ierr) 173 call CHKERR(ierr) 174 call lis_matrix_get_nnz(A,nnz,ierr) 175 176#ifdef USE_MPI 177 call MPI_Allreduce(nnz,innz,1,LIS_MPI_INTEGER,MPI_SUM,comm,ierr) 178 nnz = innz 179#endif 180 181 if ( my_rank .EQ. 0 ) then 182 write(*,'(a,i0,a,i0,a,i0,a)') 'matrix size = ', nn, ' x ', nn, ' (', nnz, ' nonzero entries)' 183 write(*,'(a)') 184 endif 185 186 call lis_matrix_duplicate(A,A0,ierr) 187 call CHKERR(ierr) 188 call lis_matrix_set_type(A0,matrix_type,ierr) 189 call lis_matrix_convert(A,A0,ierr) 190 call CHKERR(ierr) 191 call lis_matrix_destroy(A,ierr) 192 A = A0 193 194! call lis_vector_create(comm,u,ierr) 195! call lis_vector_set_size(u,0,nn,ierr) 196 197 call lis_vector_duplicate(A,u,ierr) 198 call CHKERR(ierr) 199 call lis_vector_duplicate(A,b,ierr) 200 call CHKERR(ierr) 201 call lis_vector_duplicate(A,x,ierr) 202 call CHKERR(ierr) 203 204 call lis_vector_set_all(one,u,ierr) 205 call lis_matvec(A,u,b,ierr) 206 207 call lis_solver_create(solver,ierr) 208 call CHKERR(ierr) 209 call lis_solver_set_option('-print mem',solver,ierr) 210 call lis_solver_set_optionC(solver,ierr) 211 call CHKERR(ierr) 212 213 call lis_solve(A,b,x,solver,ierr) 214 215 call CHKERR(ierr) 216 217 call lis_solver_get_iterex(solver,iter,iter_double,iter_quad,ierr) 218 call lis_solver_get_timeex(solver,time,itime,ptime,p_c_time,p_i_time,ierr) 219 call lis_solver_get_residualnorm(solver,resid,ierr) 220 call lis_solver_get_solver(solver,nsol,ierr) 221 call lis_solver_get_solvername(nsol,solvername,ierr) 222 223 if( my_rank.eq.0 ) then 224 write(*,'(a,a,i0)') solvername,': number of iterations = ', iter 225#ifndef LONG__DOUBLE 226 write(*,'(a,a,i0)') solvername,': double = ', iter_double 227 write(*,'(a,a,i0)') solvername,': quad = ', iter_quad 228#endif 229 write(*,'(a,a,e14.7e2,a)') solvername,': elapsed time = ', time, ' sec.' 230 write(*,'(a,a,e14.7e2,a)') solvername,': preconditioner = ', ptime, ' sec.' 231 write(*,'(a,a,e14.7e2,a)') solvername,': matrix creation = ', p_c_time, ' sec.' 232 write(*,'(a,a,e14.7e2,a)') solvername,': linear solver = ', itime, ' sec.' 233 write(*,'(a,a,e14.7e2)') solvername,': relative residual = ', resid 234 write(*,'(a)') '' 235 endif 236 237! write solution 238 call getarg(4,solname) 239 call lis_output_vector(x,LIS_FMT_MM,solname,ierr); 240 241! write residual history 242 call getarg(5,resname) 243 call lis_solver_output_rhistory(solver, resname, ierr) 244 245 call lis_solver_destroy(solver,ierr) 246 call lis_matrix_destroy(A,ierr) 247 call lis_vector_destroy(u,ierr) 248 call lis_vector_destroy(x,ierr) 249 call lis_vector_destroy(b,ierr) 250 251 call lis_finalize(ierr) 252 253 stop 254 end 255 256