! rbmatmodtest -- Copyright 2000 Paolo Luchini
! http://CPLcode.net/Applications/Numerical/Multigrid/
!
! test and usage example of the rbmatmod.cpl library

!( Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE. !)

! Run as, e.g.: rbmatmodtest 1 2 localhost & rbmatmodtest 2 2 localhost

USE parallel
USE rtchecks
INTEGER iproc=1, nproc=1
! iproc=number of the present node; nproc=total number of (distributed) nodes.
IF COMMANDLINE.HI>0 THEN iproc=atoi(COMMANDLINE(1)); nproc=atoi(COMMANDLINE(2))

! Open TCP sockets using the functions tcpserver and tcpclient provided in
! parallel.cpl. COMMANDLINE(3) is the hostname or ip of the previous node;
! tcpserver is listening for a similar connection from the next node.
FILE prev,next
bufsize=800
baseport=IPPORT_USERRESERVED+111
IF iproc<nproc THEN
  next=fdopen(tcpserver(baseport+iproc),"r+")
  setvbuf(next,malloc(bufsize),_IOFBF,bufsize)
END IF
IF iproc>1 THEN
  prev=fdopen(tcpclient(COMMANDLINE(3),baseport+iproc-1),"r+")
  setvbuf(prev,malloc(bufsize),_IOFBF,bufsize)
END IF
first==(prev=NULL FILE); last==(next=NULL FILE); has_terminal==last
USE rbmatmod

! Allocate some matrix and random vectors
M=10
N=10
b=2
! The unused elements of the first and last vector must be zeroed (at
! initialization only) because they will not be skipped.
REAL mat(0..N+b,-b..b), x(1..M,-b..N+b)=0, y(1..M,-b..N+b)=0, t(1..M,-b..N+b)=0
DO mat(i)=(1,-4,6,-4,1) FOR i=0 TO N
DO x(m,i)=RAND(); y(m,i)=RAND() FOR ALL m AND i=0 TO N
! make a copy for later test
savemat=mat
! parallel equivalent to LUdecomp mat (as defined in rbmat.cpl);
! modifies mat in place
LUdecompStep mat
! parallel equivalent to t=mat\x (as defined in rbmat.cpl)
DO LeftLUDivStep1[t(m),mat,x(m)] FOR ALL m
! flush output buffer before switching from writing to reading on the prev node
FlushStep1
DO LeftLUDivStep2[t(m),mat] FOR ALL m
! flush output buffer before switching from writing to reading on the next node
FlushStep2
! calculate maximum error in the solution of the linear system
REAL err=MAX ABS{x(m,i)-SUM savemat(i,j)*t(m,i+j) FOR j=-b TO b} FOR ALL m AND i=0 TO N
IF NOT first THEN err=MAX(~,BINARY REAL FROM prev)
IF NOT last THEN WRITE BINARY TO next err; FLUSH next ELSE WRITE err
! parallel scalar product between t and y
REAL sum1=SUM t(m,i)*y(m,i) FOR ALL m AND i=0 TO N
IF NOT first THEN sum1=~+BINARY REAL FROM prev
IF NOT last THEN WRITE BINARY TO next sum1; FLUSH next ELSE WRITE sum1
! parallel equivalent to y=y/mat (as defined in rbmat.cpl)
DO RightLUDivStep1[y(m),y(m),mat] FOR ALL m
FlushStep1
DO RightLUDivStep2[y(m),mat] FOR ALL m
FlushStep2
! parallel scalar product between x and y
REAL sum2=SUM x(m,i)*y(m,i) FOR ALL m AND i=0 TO N
IF NOT first THEN sum2=~+BINARY REAL FROM prev
IF NOT last THEN WRITE BINARY TO next sum2; FLUSH next ELSE WRITE sum2
! If the two products match, the test has been successful: y*mat\x = y/mat*x !
IF err<1E-7 AND ABS(1-sum2/sum1)<1E-7 THEN WRITE "success!"