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