DOUG 0.2

DOUG_utils.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 !-------------------------------------------------------
00025 module DOUG_utils
00026 
00027   use RealKind
00028   use globals
00029   use controls
00030   use mparameters
00031 
00032   implicit none
00033 
00034 #include<doug_config.h>
00035 
00036 #ifdef D_COMPLEX
00037 #define float complex
00038 #else
00039 #define float real
00040 #endif
00041 
00042   !include 'mpif.h'
00043 
00044   private :: &
00045        util_logStreamCreate,  &
00046        util_parseArgs,        &
00047        util_initMPI,          &
00048        util_finalizeMPI,      &
00049        util_printUsage,       &
00050        util_printCtrlFileInfo,&
00051        util_printVersion,     &
00052        util_actionOnCtrlArg
00053 
00054 !!$  public :: &
00055 !!$       wait_for_debugger,            &
00056 !!$       ismaster,                     &
00057 !!$       isslave,                      &
00058 !!$       DOUG_abort,                   &
00059 !!$       DOUG_quietAbort,              &
00060 !!$       DOUG_Init,                    &
00061 !!$       DOUG_Finalize,                &
00062 !!$       CtrlData_initFromFile,        &
00063 !!$       SharedCtrlData_print,         &
00064 !!$       MasterCtrlData_print,         &
00065 !!$       CtrlData_print,               &
00066 !!$       SharedCtrlData_MPItypeCreate, &
00067 !!$       SharedCtrlData_Bcast,         &
00068 !!$       length,                       &
00069 !!$       getword,                      &
00070 !!$       tolower,                      &
00071 !!$       quicksort
00072 
00073 contains
00074 
00075 
00076   subroutine doug_callme(ival,dval,indi)
00077     integer                       :: ival
00078     !real(kind=rk), intent(in out) :: dval
00079     real, intent(in out) :: dval
00080     integer, dimension(:) :: indi
00081     !integer, inten(in) :: ival
00082     !real(kind=rk) :: dval
00083     write(6,'(a,i2,a,f10.5)') '[doug_callme] : ival=',ival,&
00084         ', dval=',dval,', indi(2)=',indi
00085     ival = ival + 10
00086     dval = dval + 10.10
00087   end subroutine doug_callme
00088 
00089 
00090   !----------------------------------------
00091   ! Allows to dinamicaly attach debugger.
00092   ! Stops rank 0 and puts all the others to
00093   ! wait in a barrier.
00094   !----------------------------------------
00095   subroutine wait_for_debugger()
00096     implicit none
00097 
00098     integer :: rank, ierr
00099 
00100     call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
00101     if (rank == 0) then
00102        write(6,*) 'Waiting for debugger attachment.  Please hit enter.'
00103        pause
00104     end if
00105     call MPI_BARRIER(MPI_COMM_WORLD, ierr)
00106   end subroutine wait_for_debugger
00107 
00108 
00109   !------------------------------
00110   ! ismaster()
00111   !------------------------------
00112   function ismaster() result(res)
00113 
00114     use globals, only : D_MASTER, myrank
00115     implicit none
00116 
00117     logical :: res
00118 
00119     res = .true.
00120     if (myrank /= D_MASTER) res = .false.
00121 
00122   end function ismaster
00123 
00124 
00125   !-----------------------------
00126   ! isslave()
00127   !-----------------------------
00128   function isslave() result(res)
00129 
00130     use globals, only : D_MASTER, myrank
00131     implicit none
00132 
00133     logical :: res
00134 
00135     res = .true.
00136     if (myrank == D_MASTER) res = .false.
00137 
00138   end function isslave
00139 
00140 
00141   !-------------------------------------------
00142   ! util_logStreamCreate()
00143   !-------------------------------------------
00144   subroutine util_logStreamCreate(stream_type)
00145 
00146     use globals, only : stream, myrank, &
00147          master_stdout, slave_stdout, &
00148          D_INIT_SERIAL, D_PMASTER_LOG_FN, D_SMASTER_LOG_FN
00149 
00150     integer, intent(in), optional       :: stream_type
00151 
00152     character(150) :: fname
00153     integer        :: n, i
00154     character(150) :: ARG
00155 
00156     ! Find -q option within command line arguments.
00157     ! Specifies quiet mode when master process logs to a file
00158     ! with the name defined by D_[P/S]MASTER_LOG_FN or by parameter
00159     ! to -q option.
00160     n = iargc()
00161     do i = 1, n
00162        call getarg(i, ARG)
00163        if (trim(ARG).eq.'-q') then   ! quiet mode - master logs to a file
00164           master_stdout = .false.
00165           if (i+1 <= n) then
00166              call getarg(i+1,ARG)
00167              if ((len(trim(ARG)) /= 0).and. &
00168                   (ARG(1:1) /= '-')) then  ! a file name was specified
00169                                            ! as argument to '-q'
00170                 if (present(stream_type).and.&
00171                      (stream_type == D_INIT_SERIAL)) then
00172                    D_SMASTER_LOG_FN = trim(ARG)
00173                 else
00174                    D_PMASTER_LOG_FN = trim(ARG)
00175                 end if
00176              end if
00177              continue
00178           end if
00179        end if
00180     end do
00181 
00182     ! Open out stream(s)
00183     if (present(stream_type).and.(stream_type == D_INIT_SERIAL)) then
00184        ! In serial code log to stdout or to a file
00185        if (master_stdout) then
00186           return
00187        else
00188           write(fname,'(A)') D_SMASTER_LOG_FN
00189        end if
00190     else
00191        ! Parallel
00192        if ((ismaster().and.master_stdout).or. &
00193             (isslave().and.slave_stdout)) then
00194           return
00195        else
00196           ! open log streams to files
00197           if (myrank == 0) then
00198              write(fname,'(A)') D_PMASTER_LOG_FN
00199           else if ((myrank > 0).and.(myrank <= 9)) then
00200              write(fname,'(A4,I1)') 'log.',myrank
00201           else if ((myrank >= 10)  .and.  (myrank <= 99)) then
00202              write(fname,'(A4,I2)') 'log.',myrank
00203           else if ((myrank >= 100) .and. (myrank <= 999)) then
00204              write(fname,'(A4,I3)') 'log.',myrank
00205           end if
00206        end if
00207     end if
00208 
00209     open(stream,file=trim(fname), FORM='FORMATTED')
00210 
00211   end subroutine util_logStreamCreate
00212 
00213 
00215   subroutine util_profStreamCreate()
00216     use globals, only : pstream
00217 
00218     character(150) :: fname
00219     integer        :: n, i
00220     character(150) :: ARG
00221     logical :: profile = .FALSE.
00222 
00223     n = iargc()
00224     do i = 1, n
00225        call getarg(i, ARG)
00226        if (trim(ARG).eq.'-p') then
00227 
00228           if (i+1 <= n) then
00229              call getarg(i+1,ARG)
00230              if ((len(trim(ARG)) /= 0).and. &
00231                   (ARG(1:1) /= '-')) then  ! a file name was specified
00232                                            ! as argument to '-p'
00233                 D_PROF_FN = trim(ARG)
00234              end if
00235              continue
00236           end if
00237 
00238           write(fname,'(A,A,I2.2)') trim(D_PROF_FN), '.', myrank
00239           profile = .TRUE.
00240        end if
00241     end do
00242 
00243     if(profile) then
00244        ! open stream
00245        open(pstream,FILE=fname)
00246     else
00247        pstream = 0
00248     end if
00249   end subroutine util_profStreamCreate
00250 
00251   !----------------------------------
00252   ! DOUG_abort()
00253   !----------------------------------
00254   subroutine DOUG_abort(message, err)
00255 
00256     use globals, only: myrank, D_ERROR_STREAM, &
00257          D_INIT_TYPE, D_INIT_SERIAL
00258 
00259     character*(*)      :: message
00260     integer, optional  :: err
00261 
00262     integer            :: error = 0, ierr
00263 
00264     if (present(err)) error = err
00265 
00266     write(D_ERROR_STREAM, *)
00267     if (ismaster()) then
00268        write (D_ERROR_STREAM, 99) error, message
00269 99     format('DOUG: Fatal error: Master is bailing out with error ',i4, &
00270             /'Error message: ',a/)
00271        write(D_ERROR_STREAM, '(a)') 'doug_ended'
00272     else
00273        write (D_ERROR_STREAM,100) myrank, error, message
00274 100    format('DOUG: Fatal error: Slave ',i3,' is bailing out', &
00275             ' with error ', i4, /'Error message: ',a/)
00276        write(D_ERROR_STREAM,'(a)') 'doug_ended'
00277     endif
00278 
00279     if (D_INIT_TYPE == D_INIT_SERIAL) then
00280        stop
00281     else
00282        error = (error * 2**MPI_ABORT_ERROR_CODE_SHIFT) + 1 ! (error<<??)+1, 
00283        call MPI_ABORT(MPI_COMM_WORLD,error,ierr)
00284     end if
00285 
00286   end subroutine DOUG_abort
00287 
00288 
00289   !----------------------------------
00290   ! DOUG_quietAbort()
00291   !----------------------------------
00292   subroutine DOUG_quietAbort()
00293 
00294     use globals, only: myrank, &
00295          D_INIT_TYPE, D_INIT_SERIAL
00296 
00297     integer :: error = 0, ierr
00298 
00299     if (D_INIT_TYPE == D_INIT_SERIAL) then
00300         stop
00301      else
00302         call MPI_ABORT(MPI_COMM_WORLD,error,ierr)
00303      end if
00304 
00305    end subroutine DOUG_quietAbort
00306 
00307 
00308   !--------------------
00309   ! util_initMPI()
00310   !--------------------
00311   subroutine util_initMPI()
00312 
00313     use globals, only : numprocs, myrank, &
00314          D_MIN_PROCS_ALLOWED, D_MPI_WAS_INITED, &
00315          MPI_rkind, MPI_ckind, MPI_fkind
00316 
00317     integer       :: ierr
00318 
00319     call MPI_INITIALIZED(D_MPI_WAS_INITED, ierr) ! NB: D_MPI_WAS_INITED is of logical type
00320 
00321     if (.not.D_MPI_WAS_INITED) call MPI_INIT(ierr)
00322     call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)
00323     call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierr)
00324 
00325     ! check number of processors we are running on
00326     if (ismaster()) then
00327        if (numprocs < D_MIN_PROCS_ALLOWED) then
00328           call DOUG_abort('[util_initMPI] : Too few processors to '//&
00329                'run on. Aborting...');
00330        end if
00331     end if
00332 
00333 
00334     ! define MPI_[rcf]kind
00335     if (precision(1.0_rk) >= 15) then
00336        MPI_rkind = MPI_DOUBLE_PRECISION
00337        MPI_ckind = MPI_DOUBLE_COMPLEX
00338 #ifdef D_COMPLEX
00339        MPI_fkind = MPI_DOUBLE_COMPLEX
00340 #else
00341        MPI_fkind = MPI_DOUBLE_PRECISION
00342 #endif
00343     else
00344        MPI_rkind = MPI_REAL
00345        MPI_ckind = MPI_COMPLEX
00346 #ifdef D_COMPLEX
00347        MPI_fkind = MPI_COMPLEX
00348 #else
00349        MPI_fkind = MPI_REAL
00350 #endif
00351     end if
00352 
00353     ! define MPI_xyzkind
00354     if (precision(1.0_xyzk) >= 15) then
00355        MPI_xyzkind = MPI_DOUBLE_PRECISION
00356     else
00357        MPI_xyzkind = MPI_REAL
00358     end if
00359 
00360 
00361   end subroutine util_initMPI
00362 
00363 
00364   !------------------------------
00365   ! DOUG_Init()
00366   !------------------------------
00367   subroutine DOUG_Init(init_type)
00368 
00369     use parameters
00370     use globals, only: D_INIT_TYPE, D_INIT_PARALLEL, D_INIT_SERIAL
00371     implicit none
00372 
00373     integer, intent(in), optional :: init_type
00374 
00375     integer                       :: stream_type = D_INIT_PARALLEL
00376 
00377     call controls_init()
00378 
00379     if (present(init_type)) then
00380        ! Initialize serial DOUG (mainly for testing some drivers)
00381        if(init_type == D_INIT_SERIAL) then
00382           D_INIT_TYPE = D_INIT_SERIAL
00383           stream_type = init_type
00384        else if (init_type == D_INIT_PARALLEL) then
00385           D_INIT_TYPE = D_INIT_PARALLEL
00386           call util_initMPI()
00387        else
00388           write(6,*) 'Wrong initialization of DOUG! Aborting...'
00389           stop
00390        end if
00391     else
00392        D_INIT_TYPE = D_INIT_PARALLEL
00393        call util_initMPI()
00394     end if
00395 
00396     ! Create logging stream(s)
00397     call util_logStreamCreate(stream_type)
00398 
00399     ! Create profiling output stream
00400     call util_profStreamCreate()
00401 
00402     ! if (ismaster()) then ! MASTER
00403     !    write(stream, '(a)') 'master thread'
00404     ! else ! SLAVES
00405     !    write(stream,'(a,i3,a)') 'slave [',myrank,'] thread'
00406     ! end if
00407 
00408     ! Parse command line arguments,
00409     ! initialize DOUG controls from a file on master and
00410     ! broadcast shared controls to slaves.
00411     if ((D_INIT_TYPE == D_INIT_PARALLEL)) then
00412 
00413        if (ismaster()) then ! MASTER
00414           call util_parseArgs()
00415           write(stream, '(a)') 'master thread'
00416 
00417           call CtrlData_initFromFile()
00418           !if (D_MSGLVL > 2) &
00419           call CtrlData_print()
00420        else ! SLAVES
00421           write(stream,'(a,i3,a)') 'slave [',myrank,'] thread'
00422        end if
00423 
00424        call SharedCtrlData_MPItypeCreate()
00425        call SharedCtrlData_Bcast()
00426 
00427     end if
00428 
00429   end subroutine DOUG_Init
00430 
00431 
00432   !----------------------------
00433   ! util_finalizeMPI()
00434   !----------------------------
00435   subroutine util_finalizeMPI()
00436 
00437     use globals, only : D_MPI_WAS_INITED
00438 
00439     integer :: ierr
00440 
00441     call MPI_BARRIER(MPI_COMM_WORLD,ierr)
00442 
00443     if (.not.D_MPI_WAS_INITED) call MPI_FINALIZE(ierr)
00444 
00445   end subroutine util_finalizeMPI
00446 
00447 
00448   !-------------------------
00449   ! DOUG_Finalize()
00450   !-------------------------
00451   subroutine DOUG_Finalize()
00452 
00453     use globals, only : stream, myrank, D_INIT_TYPE, D_INIT_PARALLEL
00454     implicit none
00455 
00456     !call some_DOUG_specific_destructors()
00457 
00458     if (D_INIT_TYPE == D_INIT_PARALLEL) call util_finalizeMPI()
00459 
00460     write(stream, *)
00461     write(stream,'(A,I3,A)') 'DOUG: <',myrank,'> doug ended'
00462 
00463     close(stream)
00464 
00465   end subroutine DOUG_Finalize
00466 
00467 
00468   !-----------------------------
00469   ! Parse command line arguments
00470   !-----------------------------
00471   subroutine util_parseArgs()
00472     use globals, only: stream, D_MSGLVL, &
00473          D_CtrlFileName, master_stdout
00474 
00475     implicit none
00476 
00477     integer       :: i, n
00478     character(150) :: ARG
00479     logical       :: skipnext, CTRL_FN_SET = .false.
00480     integer       :: ierr
00481 
00482     ! write(stream,*) 'Parsing command line arguments'
00483     n = iargc()
00484     if (n == 0) then
00485        write(stream,*) 'Warning: No command line arguments specified.'
00486        write(stream,*) '         Give -h argument to see how to invoke'//&
00487             ' the program.'
00488        write(stream,*) '         Assuming DOUG control file: '//D_CtrlFileName
00489        call flush(stream)
00490        return
00491     end if
00492     ! parse arguments
00493     skipnext = .false.
00494     do i = 1, n
00495        call getarg(i, ARG)
00496        if (D_DEBUGLVL > 2) write(stream,'(a,i3,a,a)') 'argument #', i,&
00497             ' : ', trim(ARG)
00498 
00499        ! don't parse argument if it was in pair with previous one
00500        if (skipnext.eqv.(.true.)) then
00501           skipnext = .false.
00502           cycle ! skip this cycle
00503        end if
00504 
00505        if (trim(ARG).eq.'-f') then        ! control file name
00506           call getarg(i+1,ARG)
00507           skipnext = .true.
00508           if (len(trim(ARG))==0) then     ! no file name specified
00509                                           ! as argument next to '-f'
00510              write(stream,*) 'DOUG error: no control file given'
00511              call util_printUsage()
00512              call DOUG_abort('No DOUG control file specified',-1)
00513           end if
00514           D_CtrlFileName = trim(ARG)
00515           if (D_DEBUGLVL > 2) write(stream,*) &
00516                ' control file name: '//trim(D_CtrlFileName)
00517           CTRL_FN_SET = .true.
00518 
00519        else if (trim(ARG).eq.'-q') then   ! print control file parameters
00520           ! master_stdout = .false. Already set in 'util_logStreamCreate'
00521           call getarg(i+1,ARG)
00522           if ((len(trim(ARG)) /= 0).and. &
00523                   (ARG(1:1) /= '-')) then  ! a file name was specified
00524                                            ! as argument to '-q'
00525              skipnext = .true.
00526           end if
00527        else if (trim(ARG).eq.'-p') then ! Profiling, see util_profStreamCreate
00528           call getarg(i+1,ARG)
00529           if ((len(trim(ARG)) /= 0).and. &
00530                   (ARG(1:1) /= '-')) then  ! a file name was specified
00531                                            ! as argument to '-p'
00532              skipnext = .true.
00533           end if
00534        else if (trim(ARG).eq.'-i') then   ! print control file parameters
00535           call util_printCtrlFileInfo()
00536           call DOUG_quietAbort()
00537        else if (trim(ARG).eq.'-v') then   ! print version number
00538           call util_printVersion()
00539           call DOUG_quietAbort()
00540        else if (trim(ARG).eq.'-h') then   ! print usage message
00541           call util_printUsage()
00542           call DOUG_quietAbort()
00543        else if (trim(ARG).eq.'-p4pg') then   
00544           ! do not complain (for mpich ch_p4)
00545           skipnext = .true.
00546        else if (trim(ARG).eq.'-p4wd') then
00547           ! do not complain (for mpich ch_p4)
00548           skipnext = .true.
00549        else
00550           call util_printUsage()
00551           call DOUG_abort('Unrecognised command line parameter: '//&
00552                trim(ARG),-1)
00553        end if
00554 
00555     end do
00556 
00557     if ((n > 0).and.(.not.CTRL_FN_SET)) then
00558        write(stream,*) 'Warning: No -f option specified.'
00559        write(stream,*) '         Give -h argument to see how to invoke'//&
00560             ' the program.'
00561        write(stream,*) '         Assuming DOUG control file: '//D_CtrlFileName
00562        call flush(stream)
00563     end if
00564 
00565   end subroutine util_parseArgs
00566 
00567 
00568   !----------------------------------
00569   ! Print control file parameters
00570   !----------------------------------
00571   subroutine util_printCtrlFileInfo()
00572     use globals, only: stream
00573     implicit none
00574     
00575     integer :: i
00576     
00577     write(stream,*) 'List of currently supported control parameters:'
00578     do i = 1, size(ctl_words)
00579           write(stream,*) ctl_words(i)
00580     end do
00581     call flush(stream)
00582 
00583   end subroutine util_printCtrlFileInfo
00584 
00585 
00586   !-----------------------------
00587   ! Print version number
00588   !-----------------------------
00589   subroutine util_printVersion()
00590     use globals, only: stream, D_VMAJOR, D_VMINOR
00591     implicit none
00592     character(len=8) :: version
00593     character(len=1) :: str = '1'
00594 
00595     if (D_VMINOR > 9) str = '2'
00596 
00597     write(version, '(i2,a,i'//str//')') D_VMAJOR,'.',D_VMINOR
00598     write(stream,*) 'DOUG version : '//trim(version)
00599   end subroutine util_printVersion
00600 
00601 
00602   !-------------------------------
00603   ! Print message on how to use us
00604   !-------------------------------
00605   subroutine util_printUsage()
00606     use globals, only: stream
00607     implicit none
00608 
00609     write(stream,*) 'Usage: '
00610     write(stream,*) '[mpirun -np #] ./progname [[-f file.name] '//&
00611          '[-q [file.name]] [-p [file.prefix]] or [-i][-v][-h]]'
00612     write(stream,*) '                 : w/o any parameter assumes control'//&
00613          ' file is DOUG.dat'
00614     write(stream,*) ' -f file.name    : specifies DOUG control file name'
00615     write(stream,*) ' -q [file.name]  : quiet mode - master logs to'//&
00616          ' log.0/log.DOUG or to file.name'
00617     write(stream,*) ' -i              : prints out control file parameters'//&
00618          ' and exits'
00619     write(stream,*) ' -v              : prints DOUG version and exits'
00620     write(stream,*) ' -h              : prints this message and exits'
00621     write(stream,*) ' -p [file.prefix]: quiet mode - file prefix threads log to'
00622     call flush(stream)
00623   end subroutine util_printUsage
00624 
00625 
00626   !---------------------------------------------------
00627   ! Read DOUG control file and init control parameters
00628   !---------------------------------------------------
00629   subroutine CtrlData_initFromFile(CtrlFileName)
00630     use globals, only: stream, D_MSGLVL, D_CtrlFileName, sctls, mctls
00631 
00632     implicit none
00633 
00634     character*(*), optional, intent(in) :: CtrlFileName
00635 
00636     character(100)      :: ctl_fn                ! DOUG control file name
00637     integer, parameter  :: ctl_fp = 99           ! Pointer to DOUG control file
00638     character(300)      :: line, word1, word2
00639 
00640     if (present(CtrlFileName)) then
00641        ctl_fn = trim(CtrlFileName)
00642     else
00643        ctl_fn = trim(D_CtrlFileName)
00644     end if
00645 
00646     write(stream,'(a)') ' Parsing DOUG control file: '//ctl_fn
00647 
00648     open(ctl_fp, FILE=trim(ctl_fn), &
00649          STATUS='OLD', FORM='FORMATTED', ERR=999)
00650 
00651     ! Parse file:
00652 
00653     ! read stdin line by line
00654     do while(.true.)
00655 
00656        read(ctl_fp,FMT='(300a)',END=500) line
00657 
00658        ! check for comment or blank line
00659        if (len_trim(line).gt.0 .and. &
00660             ichar(line(1:1)).ne.35) then
00661 
00662           ! Break line into two pieces: control name and its argument
00663           call getword(line,1,word1)
00664           word1 = tolower(word1)
00665           call getword(line,2,word2)
00666 
00667           ! Actually fill in all control parameters.
00668           ! Action on control and its argument.
00669           call util_actionOnCtrlArg(word1, word2)
00670        end if
00671     end do ! while
00672 
00673     ! Check for not initialized critical parameters
00674     ! and initialize them to default values
00675     ! - max number of iterations:
00676     if (sctls%solve_maxiters == -1) then
00677        sctls%solve_maxiters = 100
00678     end if
00679 
00680 500 continue ! End of file reached. Close file and exit subroutine.
00681     close(ctl_fp)
00682     
00683     ! Check, if seperate RHS vector file is used.
00684         ! Give some feedback, so that this automatic checking is not so prone to user errors.
00685         inquire(file=mctls%assembled_rhs_file, exist=sctls%useAggregatedRHS)
00686         if (sctls%useAggregatedRHS) then
00687                 write(stream,*) 'Seperate aggregated RHS file is used. RHS from elemmat_rhs_file is discarded.'
00688         else
00689                 write(stream,*) 'Seperate aggregated RHS file can not be read. RHS from elemmat_rhs_file is used instead.'
00690         endif
00691     
00692     return
00693 
00694 999 call DOUG_abort('Unable to open DOUG control file: '//ctl_fn//' ', -1)
00695 
00696   end subroutine CtrlData_initFromFile
00697 
00698 
00699   !--------------------------------------------
00700   ! Action on control and its argument
00701   ! NB: D_MSGLVL and D_DEBUGLVL are reset here!
00702   !--------------------------------------------
00703   subroutine util_actionOnCtrlArg(word1, word2)
00704     use globals, only: stream, D_MSGLVL, D_DEBUGLVL, &
00705          sctls, &  ! data shared among all processes
00706          mctls     ! data belonging to master only
00707     implicit none
00708 
00709     character*(*), intent(in) :: word1, word2
00710 
00711     integer                   :: ctl_num      ! Successive number of
00712                                               ! control parameter (counter)
00713     integer                   :: i, ival
00714     logical :: lval
00715 
00716     ! find out which argument it is and set appropriatly
00717     ! WARN if it is unknown or repeated
00718     ctl_num = -1
00719     do i = 1,DCTL_NWORDS ! defined in 'controls.F90'
00720        if (trim(word1).eq.trim(ctl_words(i))) then
00721           ctl_num = i
00722        endif
00723     enddo
00724 
00725     ! begin: SHARED DATA
00726     ! solver
00727     if (ctl_num.eq.DCTL_solver) then
00728        if (sctls%solver.ne.-1) then
00729           write(6,200) trim(word1), trim(word2)
00730        else
00731           read(word2, '(i10)') sctls%solver
00732        endif
00733 
00734     ! method
00735     elseif (ctl_num.eq.DCTL_method) then
00736        if (sctls%method.ne.-1) then
00737           write(6,200) trim(word1), trim(word2)
00738        else
00739           read(word2, '(i10)') sctls%method
00740        endif
00741 
00742     ! fine method
00743     elseif (ctl_num.eq.DCTL_fine_method) then
00744        if (sctls%fine_method.ne.-1) then
00745           write(6,200) trim(word1), trim(word2)
00746        else
00747           read(word2, '(i10)') sctls%fine_method
00748        endif
00749 
00750     ! num iterations
00751     elseif (ctl_num.eq.DCTL_num_iters) then
00752        if (sctls%num_iters.ne.-1) then
00753           write(6,200) trim(word1), trim(word2)
00754        else
00755           read(word2, '(i10)') sctls%num_iters
00756        endif
00757 
00758     ! grid size
00759     elseif (ctl_num.eq.DCTL_grid_size) then
00760        if (sctls%grid_size.ne.-1) then
00761           write(6,200) trim(word1), trim(word2)
00762        else
00763           read(word2, '(i10)') sctls%grid_size
00764        endif
00765 
00766     ! coarse method
00767     elseif (ctl_num.eq.DCTL_coarse_method) then
00768        if (sctls%coarse_method.ne.-1) then
00769           write(6,200) trim(word1), trim(word2)
00770        else
00771           read(word2, '(i10)') sctls%coarse_method
00772        endif
00773 
00774     ! num_subdomains
00775     elseif (ctl_num.eq.DCTL_num_subdomains) then
00776        if (sctls%num_subdomains.ne.-1) then
00777           write(6,200) trim(word1), trim(word2)
00778        else
00779           read(word2, '(i10)') sctls%num_subdomains
00780        endif
00781 
00782     ! levels
00783     elseif (ctl_num.eq.DCTL_levels) then
00784        if (sctls%levels.ne.-1) then
00785           write(6,200) trim(word1), trim(word2)
00786        else
00787           read(word2, '(i10)') sctls%levels
00788        endif
00789 
00790     ! overlap
00791     elseif (ctl_num.eq.DCTL_overlap) then
00792        if (sctls%overlap.ne.-1) then
00793           write(6,200) trim(word1), trim(word2)
00794        else
00795           read(word2, '(i10)') sctls%overlap
00796        endif
00797 
00798     ! smoothers
00799     elseif (ctl_num.eq.DCTL_smoothers) then
00800        if (sctls%smoothers.ne.-1) then
00801           write(6,200) trim(word1), trim(word2)
00802        else
00803           read(word2, '(i10)') sctls%smoothers
00804        endif
00805 
00806     ! input_type
00807     elseif (ctl_num.eq.DCTL_input_type) then
00808        if (sctls%input_type.ne.-1) then
00809           write(6,200) trim(word1), trim(word2)
00810        else
00811           read(word2, '(i10)') sctls%input_type
00812        endif
00813 
00814     ! matrix_type
00815     elseif (ctl_num.eq.DCTL_matrix_type) then
00816        if (sctls%matrix_type.ne.-1) then
00817           write(6,200) trim(word1), trim(word2)
00818        else
00819           read(word2, '(i10)') sctls%matrix_type
00820        endif
00821 
00822     ! initial_guess
00823     elseif (ctl_num.eq.DCTL_initial_guess) then
00824        if (sctls%initial_guess.ne.-1) then
00825           write(6,200) trim(word1), trim(word2)
00826        else
00827           read(word2, '(i10)') sctls%initial_guess
00828        endif
00829 
00830     ! number_of_blocks
00831     elseif (ctl_num.eq.DCTL_number_of_blocks) then
00832        if (sctls%number_of_blocks.ne.-1) then
00833           write(6,200) trim(word1), trim(word2)
00834        else
00835           read(word2, '(i10)') sctls%number_of_blocks
00836        endif
00837 
00838     ! radius1
00839     elseif (ctl_num.eq.DCTL_radius1) then
00840        if (sctls%radius1.ne.-1) then
00841           write(6,200) trim(word1), trim(word2)
00842        else
00843           read(word2, '(i10)') sctls%radius1
00844        endif
00845 
00846     ! radius2
00847     elseif (ctl_num.eq.DCTL_radius2) then
00848        if (sctls%radius2.ne.-1) then
00849           write(6,200) trim(word1), trim(word2)
00850        else
00851           read(word2, '(i10)') sctls%radius2
00852        endif
00853 
00854     ! minasize1
00855     elseif (ctl_num.eq.DCTL_minasize1) then
00856        if (sctls%minasize1.ne.-1) then
00857           write(6,200) trim(word1), trim(word2)
00858        else
00859           read(word2, '(i10)') sctls%minasize1
00860        endif
00861 
00862     ! minasize2
00863     elseif (ctl_num.eq.DCTL_minasize2) then
00864        if (sctls%minasize2.ne.-1) then
00865           write(6,200) trim(word1), trim(word2)
00866        else
00867           read(word2, '(i10)') sctls%minasize2
00868        endif
00869 
00870     ! maxasize1
00871     elseif (ctl_num.eq.DCTL_maxasize1) then
00872        if (sctls%maxasize1.ne.-1) then
00873           write(6,200) trim(word1), trim(word2)
00874        else
00875           read(word2, '(i10)') sctls%maxasize1
00876        endif
00877 
00878     ! maxasize2
00879     elseif (ctl_num.eq.DCTL_maxasize2) then
00880        if (sctls%maxasize2.ne.-1) then
00881           write(6,200) trim(word1), trim(word2)
00882        else
00883           read(word2, '(i10)') sctls%maxasize2
00884        endif
00885 
00886     ! debug
00887     elseif (ctl_num.eq.DCTL_debug) then
00888        if (sctls%debug.ne.-1) then
00889           write(6,200) trim(word1), trim(word2)
00890        else
00891           read(word2, '(i10)') sctls%debug
00892           D_DEBUGLVL = sctls%debug
00893        endif
00894 
00895     ! verbose
00896     elseif (ctl_num.eq.DCTL_verbose) then
00897        if (sctls%verbose.ne.-1) then
00898           write(6,200) trim(word1), trim(word2)
00899        else
00900           read(word2, '(i10)') sctls%verbose
00901           D_MSGLVL = sctls%verbose
00902        endif
00903 
00904     ! plotting
00905     elseif (ctl_num.eq.DCTL_plotting) then
00906        if (sctls%plotting.ne.-1) then
00907           write(6,200) trim(word1), trim(word2)
00908        else
00909           read(word2, '(i10)') sctls%plotting
00910           ! ? = sctls%plotting
00911        endif
00912 
00913     ! strong1
00914     elseif (ctl_num.eq.DCTL_strong1) then
00915        if (sctls%strong1.ne.-1.0_rk) then
00916           write(6,200) trim(word1), trim(word2)
00917        else
00918           read(word2, '(f10.5)') sctls%strong1
00919        endif
00920 
00921     ! strong2
00922     elseif (ctl_num.eq.DCTL_strong2) then
00923        if (sctls%strong2.ne.-1.0_rk) then
00924           write(6,200) trim(word1), trim(word2)
00925        else
00926           read(word2, '(f10.5)') sctls%strong2
00927        endif
00928 
00929     ! solve_tolerance
00930     elseif (ctl_num.eq.DCTL_solve_tolerance) then
00931        if (sctls%solve_tolerance.ne.-1.0_rk) then
00932           write(6,200) trim(word1), trim(word2)
00933        else
00934           read(word2, '(f10.5)') sctls%solve_tolerance
00935        endif
00936 
00937     ! solve_maxiters
00938     elseif (ctl_num == DCTL_solve_maxiters) then
00939        if (sctls%solve_maxiters /= -1) then
00940               write(6,200) trim(word1), trim(word2)
00941        else
00942           read(word2, '(i10)') sctls%solve_maxiters
00943        endif
00944 
00945     ! symmstruct
00946     elseif (ctl_num.eq.DCTL_symmstruct) then
00947        if (sctls%symmstruct.neqv.(.false.)) then
00948           write(6,200) trim(word1), trim(word2)
00949        else
00950           read(word2, '(l1)') sctls%symmstruct
00951        endif
00952 
00953     ! symmnumeric
00954     elseif (ctl_num.eq.DCTL_symmnumeric) then
00955        if (sctls%symmnumeric.neqv.(.false.)) then
00956           write(6,200) trim(word1), trim(word2)
00957        else
00958           read(word2, '(l1)') sctls%symmnumeric
00959        endif
00960 
00961     elseif (ctl_num.eq.DCTL_interpolation_type) then
00962        if (sctls%interpolation_type.ne.-1) then
00963           write(6,200) trim(word1), trim(word2)
00964        else
00965           read(word2, '(i10)') sctls%interpolation_type
00966        endif
00967 
00968     ! end: SHARED DATA
00969 
00970     ! begin: MASTER DATA
00971     ! assembled_mtx
00972     elseif (ctl_num.eq.DCTL_assembled_mtx_file) then
00973        if (len_trim(mctls%assembled_mtx_file).ne.0) then
00974           write(6,200) trim(word1), trim(word2)
00975        else
00976           mctls%assembled_mtx_file = trim(word2)
00977        endif
00978         
00979         ! assembled_mtx_format
00980     elseif (ctl_num.eq.DCTL_assembled_mtx_format) then
00981        if (mctls%assembled_mtx_format.ne.-1) then
00982           write(6,200) trim(word1), trim(word2)
00983        else
00984           read(word2, '(i10)') mctls%assembled_mtx_format
00985        endif
00986 
00987     ! assembled_rhs
00988     elseif (ctl_num.eq.DCTL_assembled_rhs_file) then
00989        if (len_trim(mctls%assembled_rhs_file).ne.0) then
00990           write(6,200) trim(word1), trim(word2)
00991        else
00992           mctls%assembled_rhs_file = trim(word2)
00993        endif
00994         
00995         ! assembled_rhs_format
00996     elseif (ctl_num.eq.DCTL_assembled_rhs_format) then
00997        if (mctls%assembled_rhs_format.ne.-1) then
00998           write(6,200) trim(word1), trim(word2)
00999        else
01000           read(word2, '(i10)') mctls%assembled_rhs_format
01001        endif
01002 
01003     ! info_file
01004     elseif (ctl_num.eq.DCTL_info_file) then
01005        if (len_trim(mctls%info_file).ne.0) then
01006           write(6,200) trim(word1), trim(word2)
01007        else
01008           mctls%info_file = trim(word2)
01009        endif
01010 
01011     ! elemmat_rhs_file
01012     elseif (ctl_num.eq.DCTL_elemmat_rhs_file) then
01013        if (len_trim(mctls%elemmat_rhs_file).ne.0) then
01014           write(6,200) trim(word1), trim(word2)
01015        else
01016           mctls%elemmat_rhs_file = trim(word2)
01017        endif
01018 
01019     ! freedom_lists_file
01020     elseif (ctl_num.eq.DCTL_freedom_lists_file) then
01021        if (len_trim(mctls%freedom_lists_file).ne.0) then
01022           write(6,200) trim(word1), trim(word2)
01023        else
01024           mctls%freedom_lists_file = trim(word2)
01025        endif
01026 
01027     ! coords_file
01028     elseif (ctl_num.eq.DCTL_coords_file) then
01029        if (len_trim(mctls%coords_file).ne.0) then
01030           write(6,200) trim(word1), trim(word2)
01031        else
01032           mctls%coords_file = trim(word2)
01033        endif
01034 
01035     ! freemap_file
01036     elseif (ctl_num.eq.DCTL_freemap_file) then
01037        if (len_trim(mctls%freemap_file).ne.0) then
01038           write(6,200) trim(word1), trim(word2)
01039        else
01040           mctls%freemap_file = trim(word2)
01041        endif
01042 
01043     ! freedom_mask_file
01044     elseif (ctl_num.eq.DCTL_freedom_mask_file) then
01045        if (len_trim(mctls%freedom_mask_file).ne.0) then
01046           write(6,200) trim(word1), trim(word2)
01047        else
01048           mctls%freedom_mask_file = trim(word2)
01049        endif
01050 
01051     ! start_vec_file
01052     elseif (ctl_num.eq.DCTL_start_vec_file) then
01053        if (len_trim(mctls%start_vec_file).ne.0) then
01054           write(6,200) trim(word1), trim(word2)
01055        else
01056           mctls%start_vec_file = trim(word2)
01057        endif
01058 
01059     ! solution_file
01060     elseif (ctl_num.eq.DCTL_solution_file) then
01061        if (len_trim(mctls%solution_file).ne.0) then
01062           write(6,200) trim(word1), trim(word2)
01063        else
01064           mctls%solution_file = trim(word2)
01065        endif
01066 
01067     ! start_vec_type
01068     elseif (ctl_num.eq.DCTL_start_vec_type) then
01069        if (mctls%start_vec_type.ne.-1) then
01070           write(6,200) trim(word1), trim(word2)
01071        else
01072           read(word2, '(i10)') mctls%start_vec_type
01073        endif
01074 
01075     ! solution_format
01076     elseif (ctl_num.eq.DCTL_solution_format) then
01077        if (mctls%solution_format.ne.-1) then
01078           write(6,200) trim(word1), trim(word2)
01079        else
01080           read(word2, '(i10)') mctls%solution_format
01081        endif
01082     
01083     ! maxcie
01084     elseif (ctl_num.eq.DCTL_maxcie) then
01085        if (mctls%maxcie.ne.-1) then
01086           write(6,200) trim(word1), trim(word2)
01087        else
01088           read(word2, '(i10)') mctls%maxcie
01089        endif
01090 
01091     ! maxnd
01092     elseif (ctl_num.eq.DCTL_maxnd) then
01093        if (mctls%maxnd.ne.-1) then
01094           write(6,200) trim(word1), trim(word2)
01095        else
01096           read(word2, '(i10)') mctls%maxnd
01097        endif
01098 
01099     ! cutbal
01100     elseif (ctl_num.eq.DCTL_cutbal) then
01101        if (mctls%cutbal.ne.-1) then
01102           write(6,200) trim(word1), trim(word2)
01103        else
01104           read(word2, '(i10)') mctls%cutbal
01105        endif
01106 
01107     ! center_type
01108     elseif (ctl_num.eq.DCTL_center_type) then
01109        if (mctls%center_type.ne.-1) then
01110           write(6,200) trim(word1), trim(word2)
01111        else
01112           read(word2, '(i10)') mctls%center_type
01113        endif
01114 
01115     ! hanging_nodes
01116     elseif (ctl_num.eq.DCTL_hanging_nodes) then
01117        if (mctls%hanging_nodes.neqv.(.false.)) then
01118           write(6,200) trim(word1), trim(word2)
01119        else
01120           read(word2, '(l1)') mctls%hanging_nodes
01121        endif
01122 
01123     ! dump_matrix_only
01124     elseif (ctl_num.eq.DCTL_dump_matrix_only) then
01125        if (mctls%dump_matrix_only.neqv.(.false.)) then
01126           write(6,200) trim(word1), trim(word2)
01127        else
01128           read(word2, '(l1)') mctls%dump_matrix_only
01129        endif
01130 
01131         ! dump_matrix_file
01132     elseif (ctl_num.eq.DCTL_dump_matrix_file) then
01133        if (len_trim(mctls%dump_matrix_file).ne.0) then
01134           write(6,200) trim(word1), trim(word2)
01135        else
01136           mctls%dump_matrix_file = trim(word2)
01137        endif
01138        
01139     ! end: MASTER DATA
01140 
01141     ! nothing was recognised
01142     else
01143        write(6,210) trim(word1), trim(word2)
01144     end if
01145 
01146 200 format('Warning: Multiple use of: [',a,'] argument: ',a)
01147 210 format('Warning: Unrecognised control: [',a,'] argument: ',a)
01148 
01149   end subroutine util_actionOnCtrlArg
01150 
01151 
01152   !---------------------------------------------------------------------
01153   ! Print data type representing control data shared among all processes
01154   !---------------------------------------------------------------------
01155   subroutine SharedCtrlData_print(noheader)
01156     use globals, only: stream, &
01157          sctls ! data shared among all processes
01158     implicit none
01159 
01160     logical, optional, intent(in) :: noheader
01161 
01162     integer       :: j, lens
01163     character(24) :: fmta, fmti, fmtl, fmtr, fmtc
01164 
01165     lens = length(ctl_words(1))
01166     do j = 2,DCTL_NWORDS
01167        if (lens < length(ctl_words(j))) then
01168           lens = length(ctl_words(j))
01169        end if
01170     end do
01171 
01172     write(fmta, fmt="('(a',i5,','' = ')") lens
01173     fmtl = fmta(1:length(fmta))//" ',l5)"
01174     fmti = fmta(1:length(fmta))//" ',i5)"
01175     fmtr = fmta(1:length(fmta))//" ',1p,e12.3)"
01176     fmtc = fmta(1:length(fmta))//" ',a)"
01177 
01178     if (.not.present(noheader)) then
01179        write(stream,*)
01180        write(stream,*) 'Shared control parameters:'
01181     end if
01182     write(stream,fmti) &
01183          ctl_words(DCTL_solver)(1:length(ctl_words(DCTL_solver))), &
01184          sctls%solver
01185     write(stream,fmti) &
01186          ctl_words(DCTL_method)(1:length(ctl_words(DCTL_method))), &
01187          sctls%method
01188     write(stream,fmti) &
01189          ctl_words(DCTL_fine_method)(1:length(ctl_words(DCTL_fine_method))), &
01190          sctls%fine_method
01191     write(stream,fmti) &
01192          ctl_words(DCTL_num_iters)(1:length(ctl_words(DCTL_num_iters))), &
01193          sctls%num_iters
01194     write(stream,fmti) &
01195          ctl_words(DCTL_grid_size)(1:length(ctl_words(DCTL_grid_size))), &
01196          sctls%grid_size
01197     write(stream,fmti) &
01198          ctl_words(DCTL_coarse_method)(1:length(ctl_words(DCTL_coarse_method))), &
01199          sctls%coarse_method
01200     write(stream,fmti) &
01201          ctl_words(DCTL_num_subdomains)(1:length(ctl_words(DCTL_num_subdomains))), &
01202          sctls%num_subdomains
01203     write(stream,fmti) &
01204          ctl_words(DCTL_levels)(1:length(ctl_words(DCTL_levels))), &
01205          sctls%levels
01206     write(stream,fmti) &
01207          ctl_words(DCTL_overlap)(1:length(ctl_words(DCTL_overlap))), &
01208          sctls%overlap
01209     write(stream,fmti) &
01210          ctl_words(DCTL_smoothers)(1:length(ctl_words(DCTL_smoothers))), &
01211          sctls%smoothers
01212     write(stream,fmtl) &
01213          ctl_words(DCTL_symmstruct)(1:length(ctl_words(DCTL_symmstruct))), &
01214          sctls%symmstruct
01215     write(stream,fmtl) &
01216          ctl_words(DCTL_symmnumeric)(1:length(ctl_words(DCTL_symmnumeric))), &
01217          sctls%symmnumeric
01218     write(stream,fmti) &
01219          ctl_words(DCTL_input_type)(1:length(ctl_words(DCTL_input_type))), &
01220          sctls%input_type
01221     write(stream,fmti) &
01222          ctl_words(DCTL_matrix_type)(1:length(ctl_words(DCTL_matrix_type))), &
01223          sctls%matrix_type
01224     write(stream,fmti) &
01225          ctl_words(DCTL_number_of_blocks) &
01226          (1:length(ctl_words(DCTL_number_of_blocks))), sctls%number_of_blocks
01227     write(stream,fmtr) &
01228          ctl_words(DCTL_strong1) &
01229          (1:length(ctl_words(DCTL_strong1))), sctls%strong1
01230     write(stream,fmtr) &
01231          ctl_words(DCTL_strong2) &
01232          (1:length(ctl_words(DCTL_strong2))), sctls%strong2
01233     write(stream,fmtr) &
01234          ctl_words(DCTL_solve_tolerance) &
01235          (1:length(ctl_words(DCTL_solve_tolerance))), sctls%solve_tolerance
01236     write(stream,fmti) &
01237          ctl_words(DCTL_solve_maxiters) &
01238          (1:length(ctl_words(DCTL_solve_maxiters))), sctls%solve_maxiters
01239     write(stream,fmti) &
01240          ctl_words(DCTL_radius1)(1:length(ctl_words(DCTL_radius1))), &
01241          sctls%radius1
01242     write(stream,fmti) &
01243          ctl_words(DCTL_radius2)(1:length(ctl_words(DCTL_radius2))), &
01244          sctls%radius2
01245     write(stream,fmti) &
01246          ctl_words(DCTL_minasize1)(1:length(ctl_words(DCTL_minasize1))), &
01247          sctls%minasize1
01248     write(stream,fmti) &
01249          ctl_words(DCTL_minasize2)(1:length(ctl_words(DCTL_minasize2))), &
01250          sctls%minasize2
01251     write(stream,fmti) &
01252          ctl_words(DCTL_maxasize1)(1:length(ctl_words(DCTL_maxasize1))), &
01253          sctls%maxasize1
01254     write(stream,fmti) &
01255          ctl_words(DCTL_maxasize2)(1:length(ctl_words(DCTL_maxasize2))), &
01256          sctls%maxasize2
01257     write(stream,fmti) &
01258          ctl_words(DCTL_debug)(1:length(ctl_words(DCTL_debug))), &
01259          sctls%debug
01260     write(stream,fmti) &
01261          ctl_words(DCTL_verbose)(1:length(ctl_words(DCTL_verbose))), &
01262          sctls%verbose
01263     write(stream,fmti) &
01264          ctl_words(DCTL_plotting)(1:length(ctl_words(DCTL_plotting))), &
01265          sctls%plotting
01266     write(stream,fmti) &
01267          ctl_words(DCTL_initial_guess) &
01268          (1:length(ctl_words(DCTL_initial_guess))), sctls%initial_guess
01269     write(stream,fmti) &
01270          ctl_words(DCTL_interpolation_type)&
01271              (1:length(ctl_words(DCTL_interpolation_type))), &
01272          sctls%interpolation_type
01273     call flush(stream)
01274 
01275   end subroutine SharedCtrlData_print
01276 
01277 
01278   !-------------------------------------------------------------------
01279   ! Print data type representing control data belonging to master only
01280   !-------------------------------------------------------------------
01281   subroutine MasterCtrlData_print(noheader)
01282    use globals, only: stream, &
01283          mctls ! control data known to master only
01284     implicit none
01285 
01286     logical, optional, intent(in) :: noheader
01287 
01288     integer       :: j, lens
01289     character(24) :: fmta, fmti, fmtr, fmtc, fmtl
01290 
01291     lens = length(ctl_words(1))
01292     do j = 2,DCTL_NWORDS
01293        if (lens < length(ctl_words(j))) then
01294           lens = length(ctl_words(j))
01295        end if
01296     end do
01297 
01298     write(fmta, fmt="('(a',i5,','' = ')") lens
01299     fmti = fmta(1:length(fmta))//" ',i5)"
01300     fmtr = fmta(1:length(fmta))//" ',1p,e12.3)"
01301     fmtc = fmta(1:length(fmta))//" ',a)"
01302     fmtl = fmta(1:length(fmta))//" ',l1)"
01303 
01304     if (.not.present(noheader)) then
01305        write(stream,*)
01306        write(stream,*) 'Master control parameters:'
01307     end if
01308     write(stream,fmti) &
01309          ctl_words(DCTL_solution_format) &
01310          (1:length(ctl_words(DCTL_solution_format))), mctls%solution_format
01311     write(stream,fmti) &
01312          ctl_words(DCTL_start_vec_type) &
01313          (1:length(ctl_words(DCTL_start_vec_type))), mctls%start_vec_type
01314     write(stream,fmti) &
01315          ctl_words(DCTL_assembled_rhs_format) &
01316          (1:length(ctl_words(DCTL_assembled_rhs_format))), & 
01317          mctls%assembled_rhs_format
01318     write(stream,fmtc) &
01319          ctl_words(DCTL_assembled_rhs_file) &
01320          (1:length(ctl_words(DCTL_assembled_rhs_file))), &
01321          trim(mctls%assembled_rhs_file)
01322     if (sctls%input_type==DISTRIBUTION_TYPE_ASSEMBLED) then
01323       write(stream,fmtc) &
01324            ctl_words(DCTL_assembled_mtx_file) &
01325            (1:length(ctl_words(DCTL_assembled_mtx_file))), &
01326            trim(mctls%assembled_mtx_file)
01327       write(stream,fmti) &
01328            ctl_words(DCTL_assembled_mtx_format) &
01329            (1:length(ctl_words(DCTL_assembled_mtx_format))), &
01330            mctls%assembled_mtx_format
01331     elseif (sctls%input_type==DISTRIBUTION_TYPE_ELEMENTAL) then
01332       write(stream,fmtc) &
01333            ctl_words(DCTL_info_file)(1:length(ctl_words(DCTL_info_file))), &
01334            trim(mctls%info_file)
01335       write(stream,fmtc) &
01336            ctl_words(DCTL_freedom_lists_file) &
01337            (1:length(ctl_words(DCTL_freedom_lists_file))), &
01338            trim(mctls%freedom_lists_file)
01339       write(stream,fmtc) &
01340            ctl_words(DCTL_elemmat_rhs_file) &
01341            (1:length(ctl_words(DCTL_elemmat_rhs_file))), &
01342            trim(mctls%elemmat_rhs_file)
01343       write(stream,fmtc) &
01344            ctl_words(DCTL_coords_file)(1:length(ctl_words(DCTL_coords_file))),&
01345            trim(mctls%coords_file)
01346       write(stream,fmtc) &
01347            ctl_words(DCTL_freemap_file) &
01348            (1:length(ctl_words(DCTL_freemap_file))), trim(mctls%freemap_file)
01349       write(stream,fmtc) &
01350            ctl_words(DCTL_freedom_mask_file)&
01351            (1:length(ctl_words(DCTL_freedom_mask_file))), &
01352            trim(mctls%freedom_mask_file)
01353       write(stream,fmtc) &
01354            ctl_words(DCTL_solution_file) &
01355            (1:length(ctl_words(DCTL_solution_file))), trim(mctls%solution_file)
01356       write(stream,fmtc) &
01357            ctl_words(DCTL_start_vec_file) &
01358            (1:length(ctl_words(DCTL_start_vec_file))), &
01359            trim(mctls%start_vec_file)
01360       write(stream,fmtc) &
01361            ctl_words(DCTL_dump_matrix_file) &
01362            (1:length(ctl_words(DCTL_dump_matrix_file))), &
01363            trim(mctls%dump_matrix_file)
01364       write(stream,fmtl) &
01365            ctl_words(DCTL_dump_matrix_only) &
01366            (1:length(ctl_words(DCTL_dump_matrix_only))), mctls%dump_matrix_only
01367 
01368       write(stream,fmti) &
01369          ctl_words(DCTL_maxcie) &
01370          (1:length(ctl_words(DCTL_maxcie))), mctls%maxcie
01371       write(stream,fmti) &
01372          ctl_words(DCTL_maxnd) &
01373          (1:length(ctl_words(DCTL_maxnd))), mctls%maxnd
01374       write(stream,fmti) &
01375          ctl_words(DCTL_cutbal) &
01376          (1:length(ctl_words(DCTL_cutbal))), mctls%cutbal
01377 
01378 
01379 
01380 
01381     endif
01382     call flush(stream)
01383 
01384   end subroutine MasterCtrlData_print
01385 
01386 
01387   !---------------------------------------------------
01388   ! Print out DOUG control parameters and their values
01389   !---------------------------------------------------
01390   subroutine CtrlData_print()
01391     use globals, only: stream
01392     implicit none
01393 
01394     logical :: noheader=.true.
01395 
01396     write(stream,*) 'Control parameters:'
01397     call SharedCtrlData_print(noheader)
01398     call MasterCtrlData_print(noheader)
01399 
01400   end subroutine CtrlData_print
01401 
01402 
01403   !-------------------------------------------------------
01404   ! Create MPI type for 'SharedCtrlData' derived data type
01405   !-------------------------------------------------------
01406   subroutine SharedCtrlData_MPItypeCreate()
01407     use globals, only: sctls, D_MPI_SCTLS_TYPE, MPI_rkind
01408     implicit none
01409 
01410     integer, parameter          :: nblocks = 31  ! Number of type components
01411     integer, dimension(nblocks) :: types,        
01412                                    blocklengths, 
01413                                    addresses,    
01414                                    displacements
01415     integer                     :: i, ierr
01416 
01417 
01418     types = (/MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, &
01419          MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, &
01420          MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, &
01421          MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, &
01422          MPI_INTEGER, &
01423          MPI_INTEGER, MPI_rkind, MPI_rkind, MPI_rkind, &
01424          MPI_INTEGER, MPI_LOGICAL, &
01425          MPI_LOGICAL, &
01426          MPI_INTEGER, MPI_LOGICAL, &
01427          MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, &
01428          MPI_INTEGER /)
01429 
01430     blocklengths = (/1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
01431                      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1/)
01432 
01433     call MPI_ADDRESS(sctls%solver,           addresses( 1), ierr)
01434     call MPI_ADDRESS(sctls%method,           addresses( 2), ierr)
01435     call MPI_ADDRESS(sctls%levels,           addresses( 3), ierr)
01436     call MPI_ADDRESS(sctls%overlap,          addresses( 4), ierr)
01437     call MPI_ADDRESS(sctls%smoothers,        addresses( 5), ierr)
01438     call MPI_ADDRESS(sctls%input_type,       addresses( 6), ierr)
01439     call MPI_ADDRESS(sctls%matrix_type,      addresses( 7), ierr)
01440     call MPI_ADDRESS(sctls%initial_guess,    addresses( 8), ierr)
01441     call MPI_ADDRESS(sctls%number_of_blocks, addresses( 9), ierr)
01442     call MPI_ADDRESS(sctls%radius1,          addresses(10), ierr)
01443     call MPI_ADDRESS(sctls%radius2,          addresses(11), ierr)
01444     call MPI_ADDRESS(sctls%minasize1,        addresses(12), ierr)
01445     call MPI_ADDRESS(sctls%minasize2,        addresses(13), ierr)
01446     call MPI_ADDRESS(sctls%maxasize1,        addresses(14), ierr)
01447     call MPI_ADDRESS(sctls%maxasize2,        addresses(15), ierr)
01448     call MPI_ADDRESS(sctls%debug,            addresses(16), ierr)
01449     call MPI_ADDRESS(sctls%verbose,          addresses(17), ierr)
01450     call MPI_ADDRESS(sctls%plotting,         addresses(18), ierr)
01451     call MPI_ADDRESS(sctls%strong1,          addresses(19), ierr)
01452     call MPI_ADDRESS(sctls%strong2,          addresses(20), ierr)
01453     call MPI_ADDRESS(sctls%solve_tolerance,  addresses(21), ierr)
01454     call MPI_ADDRESS(sctls%solve_maxiters,   addresses(22), ierr)
01455     call MPI_ADDRESS(sctls%symmstruct,       addresses(23), ierr)
01456     call MPI_ADDRESS(sctls%symmnumeric,      addresses(24), ierr)
01457     call MPI_ADDRESS(sctls%interpolation_type, addresses(25), ierr)
01458     call MPI_ADDRESS(sctls%useAggregatedRHS, addresses(26), ierr)
01459     call MPI_ADDRESS(sctls%coarse_method,    addresses(27), ierr)
01460     call MPI_ADDRESS(sctls%num_subdomains,    addresses(28), ierr)
01461     call MPI_ADDRESS(sctls%fine_method,    addresses(29), ierr)
01462     call MPI_ADDRESS(sctls%num_iters,    addresses(30), ierr)
01463     call MPI_ADDRESS(sctls%grid_size,    addresses(31), ierr)
01464 
01465     do i = 1,nblocks
01466        displacements(i) = addresses(i) - addresses(1)
01467     end do
01468 
01469     call MPI_TYPE_STRUCT(nblocks, blocklengths, displacements, types, &
01470          D_MPI_SCTLS_TYPE, ierr);
01471     call MPI_TYPE_COMMIT(D_MPI_SCTLS_TYPE, ierr)
01472 
01473   end subroutine SharedCtrlData_MPItypeCreate
01474 
01475 
01476   !------------------------------------------------------
01477   ! Broadcast shared control data
01478   ! NB: D_MSGLVL and D_DEBUGLVL are reset here on slaves!
01479   !------------------------------------------------------
01480   subroutine SharedCtrlData_Bcast()
01481     implicit none
01482 
01483     integer :: ierr
01484 
01485     call MPI_BCAST(sctls, 1, D_MPI_SCTLS_TYPE,&
01486          D_MASTER, MPI_COMM_WORLD, ierr)
01487 
01488     ! Initialise messaging and debug levels on slaves
01489     if (isslave()) then
01490        D_MSGLVL   = sctls%verbose
01491        D_DEBUGLVL = sctls%debug
01492     end if
01493 
01494   end subroutine SharedCtrlData_Bcast
01495 
01496 
01497   !------------------------------------------------
01498   ! Length of a string without back trailing spaces
01499   !------------------------------------------------
01500   function length(string)
01501     implicit none
01502     character*(*), intent(in) :: string
01503     integer                   :: length, i
01504 
01505     do i = len(string),1,-1
01506        if (string(i:i).ne.' ') then
01507           length = i
01508           return
01509        endif
01510     enddo
01511     length = 0
01512   end function length
01513 
01514 
01515   !-------------------------------------
01516   ! Find the num'th word in a string
01517   !-------------------------------------
01518   subroutine getword(string, num, wordx)
01519     implicit none
01520 
01521     character*(*), intent(in)  :: string
01522     character*(*), intent(out) :: wordx
01523     integer                    :: num,i,j,k,n
01524 
01525     n = len_trim(string)
01526 
01527     j = 1
01528     do i = 1,num
01529        do while (string(j:j).eq.' '.and.j.le.n)
01530           j = j+1
01531        enddo
01532        if (j.gt.n) then
01533           wordx=' '
01534           return
01535        endif
01536 
01537        k = j
01538        do while (string(j:j).ne.' ')
01539           j = j+1
01540        enddo
01541     enddo
01542 
01543     wordx = string(k:j)
01544   end subroutine getword
01545 
01546 
01547   !-------------------------------------
01548   ! convert a string to lower case
01549   !-------------------------------------
01550   function tolower(string) result(lower)
01551     implicit none
01552 
01553     character*(*), intent(in) :: string
01554     character*(300)           :: lower
01555     integer                   :: i
01556 
01557     do i=1,len_trim(string)
01558        if (ichar(string(i:i)).ge.65 .and. &
01559             ichar(string(i:i)).le.90) then
01560           lower(i:i)=char(ichar(string(i:i))+32)
01561        else
01562           lower(i:i)=string(i:i)
01563        endif
01564     enddo
01565 
01566     do i=len_trim(string)+1,len_trim(lower)
01567        lower(i:i)=' '
01568     enddo
01569 
01570   end function tolower
01571   
01572   !------------------------------------------------
01575   subroutine FindFreeIOUnit(found,iounit)
01576         implicit none
01577         
01578         logical, intent(out)  :: found        !< has a free IO unit been found?
01579         integer, intent(out)  :: iounit       !< number of free IO unit found
01580         logical               :: exi, ope
01581         found = .true.
01582         do iounit=7,99
01583            inquire(unit=iounit,exist=exi) 
01584            if( exi ) inquire(unit=iounit,opened=ope)
01585            if( exi.and..not.ope) &
01586               exit
01587            if( iounit .eq. 99) then
01588               write(stream, *) 'Could not find free I/O-Unit!'
01589               found = .false.
01590            endif
01591         enddo
01592         
01593   end subroutine FindFreeIOUnit
01594 
01595   !-------------------------------------
01598   subroutine WriteSolutionToFile(x) !TODO: Refactor, XDR stuff finds free IO unit by itself.
01599                                                                         !          open should be made inside specific subroutines.
01600         implicit none
01601            
01602     float(kind=rk), dimension(:), intent(in) :: x        !< vector to write out
01603  
01604     integer :: iounit,opened
01605         logical :: found
01606 
01607         call FindFreeIOUnit(found, iounit)
01608         if (found) then
01609    
01610       ! Text format.   
01611       if (mctls%solution_format == 0) then
01612         open(unit=iounit,iostat=opened,file=mctls%solution_file)
01613         if (opened.eq.0) &
01614           call WriteSolutionTextualFormat(iounit, x)
01615       
01616       ! Binary format.
01617       elseif (mctls%solution_format == 1) then
01618             open(unit=iounit,iostat=opened,form='unformatted',file=mctls%solution_file)
01619             if (opened.eq.0) &
01620           call WriteSolutionBinaryFormat(iounit, x)
01621 
01622 #ifdef HAVE_LIBFXDR
01623       ! XDR format 
01624       elseif (mctls%solution_format == 2) then
01625         call WriteSolutionXDRFormat(mctls%solution_file, x)
01626 #endif
01627         
01628       ! Default case.
01629       else
01630         write (stream,*) 'Solution format is not valid. Set it in the control file.'
01631         write (stream,*) 'Solution not written.'
01632       endif
01633       
01634       close(iounit)
01635     endif
01636     
01637   end subroutine WriteSolutionToFile
01638 
01639   !----------------------------------
01643   subroutine WriteSolutionTextualFormat(iounit, x)
01644         implicit none
01645            
01646         integer, intent(in)                      :: iounit   !< IO-unit to write to
01647     float(kind=rk), dimension(:), intent(in) :: x        !< vector to write out
01648     integer :: i,n
01649     
01650     n = size(x)
01651     write(iounit,'(i6)') n
01652     do i = 1,n
01653        write(iounit, '(e21.14)') x(i)
01654     end do
01655  
01656   end subroutine
01657   
01658   !----------------------------------
01666   subroutine WriteSolutionBinaryFormat(iounit, x)
01667         implicit none
01668            
01669         integer, intent(in)                      :: iounit   !< IO-unit to write to
01670     float(kind=rk), dimension(:), intent(in) :: x        !< vector to write out
01671     integer :: k
01672 
01673     print *, iounit, size(x)
01674     write(iounit) (x(k), k = 1,size(x))
01675 
01676   end subroutine
01677 
01678 #ifdef HAVE_LIBFXDR
01679   !----------------------------------
01687   subroutine WriteSolutionXDRFormat(fname, x)
01688         implicit none
01689          include 'fxdr.inc'
01690         
01691         character*(*),                intent(in) :: fname    !< file name   
01692     float(kind=rk), dimension(:), intent(in) :: x        !< vector to write to file
01693     integer :: k, fhandle, status
01694 
01695     fhandle = initxdr(fname, 'w', .FALSE.)
01696     !status = ixdrint( fhandle, size(x) )
01697         status = ixdrdmat( fhandle, size(x), x )
01698         status = ixdrclose( fhandle )
01699 
01700   end subroutine
01701 
01702 #endif
01703 
01704   !-------------------------------------
01707   subroutine quicksort(n,indx)
01708       implicit none
01709       integer :: n !< size of the array indx
01710       integer :: indx(n) !< values to be sorted
01711       !------------------------
01712       integer,parameter :: M=7,NSTACK=64
01713       integer :: i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
01714 
01715       jstack=0
01716       l=1
01717       ir=n
01718  1    if(ir-l.lt.M) then
01719          do j=l+1,ir
01720             indxt=indx(j)
01721             do i=j-1,l,-1
01722                if(indx(i).le.indxt) goto 2
01723                indx(i+1)=indx(i)
01724             enddo
01725             i=l-1
01726  2          indx(i+1)=indxt
01727          enddo
01728          if(jstack.eq.0) return
01729         ir=istack(jstack)
01730          l=istack(jstack-1)
01731          jstack=jstack-2
01732       else
01733          k=(l+ir)/2
01734          itemp=indx(k)
01735          indx(k)=indx(l+1)
01736          indx(l+1)=itemp
01737          if(indx(l).gt.indx(ir)) then
01738             itemp=indx(l)
01739             indx(l)=indx(ir)
01740             indx(ir)=itemp
01741          endif
01742          if(indx(l+1).gt.indx(ir)) then
01743             itemp=indx(l+1)
01744             indx(l+1)=indx(ir)
01745             indx(ir)=itemp
01746          endif
01747         if(indx(l).gt.indx(l+1)) then
01748             itemp=indx(l)
01749             indx(l)=indx(l+1)
01750             indx(l+1)=itemp
01751          endif
01752          i=l+1
01753          j=ir
01754          indxt=indx(l+1)
01755  3       continue
01756          i=i+1
01757          if(indx(i).lt.indxt) goto 3
01758  4       continue
01759          j=j-1
01760          if(indx(j).gt.indxt) goto 4
01761          if(j.lt.i) goto 5
01762          itemp=indx(i)
01763          indx(i)=indx(j)
01764          indx(j)=itemp
01765          goto 3
01766 5       indx(l+1)=indx(j)
01767          indx(j)=indxt
01768          jstack=jstack+2
01769          if(jstack.gt.NSTACK) then
01770             write(stream,200) NSTACK
01771  200        format('Quicksort: NSTACK=',i4,' apparently ',&
01772                  'too small for this problem')
01773             call DOUG_abort('Quicksort failed',50)
01774          endif
01775          if(ir-i+1.ge.j-l) then
01776             istack(jstack)=ir
01777             istack(jstack-1)=i
01778             ir=j-1
01779          else
01780             istack(jstack)=j-1
01781             istack(jstack-1)=l
01782             l=i
01783          endif
01784       endif
01785       goto 1
01786   end subroutine quicksort
01787 
01788 
01789   !-------------------------------------
01792   function random_integer(mini,maxi) result(irand)
01793     integer,intent(in) :: mini,maxi
01794     real :: r
01795     integer :: irand
01796     ! call random_seed before first call !
01797     call random_number(r)
01798     r=r*(maxi-mini+1)
01799     irand=mini+floor(r)
01800   end function random_integer
01801 
01802   !-------------------------------------
01805   subroutine rand_intarr(n,riarr,mini,maxi)
01806     integer,intent(in) :: n,mini,maxi
01807     integer,dimension(n),intent(out) :: riarr
01808     real,dimension(:),allocatable :: ra
01809     ! call random_seed before first call !
01810     allocate(ra(n))
01811     call random_number(ra)
01812     ra=ra*(maxi-mini+1)
01813     riarr=mini+floor(ra)
01814     deallocate(ra)
01815   end subroutine rand_intarr
01816 
01817   !-------------------------------------
01820   subroutine random_permutation(n,inds)
01821     integer,intent(in) :: n
01822     integer,dimension(:),pointer :: inds
01823     real,dimension(:),allocatable :: ra
01824     real :: rrand
01825     integer :: i,irand,k
01826     ! call random_seed before first call !
01827     call random_seed
01828     allocate(ra(n-1))
01829     call random_number(ra)
01830     ! Loop over the indeces
01831     do i=1,n-1
01832       ! take a random number in [i,n]:
01833       rrand=ra(i)*(n-i+1)
01834       irand=i+floor(rrand)
01835       k=inds(i)
01836       inds(i)=inds(irand)
01837       inds(irand)=k
01838     enddo
01839     deallocate(ra)
01840   end subroutine random_permutation
01841 
01842 end module DOUG_utils