text:
!*******************************************************************************
!>  @file bmw.f
!>  @brief Contains the main routines for Biot-Savart Magnetic Vmec
!>  Vector-Potential (BMW) code.
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!>  BMW is a code for extending fields belond the VMEC domain in a manner that
!>  ensures divergence free fields. BMW does a Biot-Savart volume integration of
!>  of the equilibrium current density to obtain a continous vector potential
!>  every where on the mgrid grid.
!>
!>  @author Mark Cianciosa
!>
!>  Below is a brief discription of the major top level objects of the code. For
!>  discriptions of lower level objects consult the referenced top level
!>  objects.
!>
!>  @ref m_grid        Object containing the vaccum field information.
!>  @ref primed_grid   Object containing the plasma currents and primed grid
!>                     positions.
!>  @ref unprimed_grid Object containing the plasma vector potential response.
!>  @ref bmw_context   Defines the main pointers to all allocated memory and
!>                     objects.
!*******************************************************************************
!*******************************************************************************
!  MAIN PROGRAM
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief BMW main program.
!>
!>  Highest level BMW routine. This computes the vector potential on the mgrid
!>  grid.
!-------------------------------------------------------------------------------
      PROGRAM bmw
      USE bmw_context
      USE bmw_commandline_parser
      USE, INTRINSIC :: iso_fortran_env, Only : output_unit
      USE bmw_state_flags
      USE safe_open_mod

      IMPLICIT NONE

!  local variables
      CLASS (bmw_parallel_context_class), POINTER :: parallel => null()
      CLASS (bmw_context_class), POINTER            :: context => null()
      CLASS (bmw_commandline_parser_class), POINTER ::                         &
     &   cl_parser => null()
      INTEGER                                       :: flags
      INTEGER                                       :: num_p
      INTEGER                                       :: io_unit
      INTEGER                                       :: status
      REAL (rprec)                                  :: start_time

!  Start of executable code
#if defined (MPI_OPT)
      CALL MPI_INIT(status)
#endif
      CALL profiler_construct

      start_time = profiler_get_start_time()

      parallel => bmw_parallel_context_class(                                  &
#if defined (MPI_OPT)
     &                                       MPI_COMM_WORLD                    &
#endif
     &                                       )

      cl_parser => bmw_commandline_parser_class(parallel)

!  Check if the required flags are set.
!      IF (.not.cl_parser%is_flag_set('-mgridf')) THEN
!         WRITE (*,1001) '-mgridf'
!         CALL bmw_commandline_parser_print_help
!      END IF
      IF (.not.cl_parser%is_flag_set('-woutf')) THEN
         WRITE (*,1001) '-woutf'
         CALL bmw_commandline_parser_print_help
      END IF
      IF (.not.cl_parser%is_flag_set('-outf')) THEN
         WRITE (*,1001) '-outf'
         CALL bmw_commandline_parser_print_help
      END IF

      CALL parallel%set_threads(cl_parser%get('-para', 1))

      flags = bmw_state_flags_off
      IF (cl_parser%is_flag_set('-force')) THEN
         flags = IBSET(flags, bmw_state_flags_force)
      END IF
      IF (cl_parser%is_flag_set('-ju')) THEN
         flags = IBSET(flags, bmw_state_flags_ju)
      END IF
      IF (cl_parser%is_flag_set('-jv')) THEN
         flags = IBSET(flags, bmw_state_flags_jv)
         flags = IBCLR(flags, bmw_state_flags_ju)
      END IF
      IF (cl_parser%is_flag_set('-siestaf')) THEN
         flags = IBSET(flags, bmw_state_flags_siesta)
         flags = IBCLR(flags, bmw_state_flags_ju)
         flags = IBCLR(flags, bmw_state_flags_jv)
      END IF

      flags = IBSET(flags, bmw_state_flags_mgrid)
      num_p = 0

      io_unit = output_unit
      IF (cl_parser%is_flag_set('-logf')) THEN
         CALL safe_open(io_unit, status, cl_parser%get_string('-logf'),        &
     &                  'replace', 'formatted', delim_in='none')
      END IF

      IF (parallel%offset .eq. 0) THEN
         WRITE (io_unit,*)
         WRITE (io_unit,1000) series
         WRITE (io_unit,*)
      END IF
      CALL parallel%report(io_unit)

      IF (cl_parser%is_flag_set('-mgridf')) THEN
         context => bmw_context_class(cl_parser%get('-mgridf'),                &
     &                                cl_parser%get('-woutf'),                 &
     &                                cl_parser%get('-siestaf'),               &
     &                                flags, num_p, parallel, io_unit)
      ELSE
         num_p = cl_parser%get('-num_p', 1)
         context => bmw_context_class(cl_parser%get('-woutf'),                 &
     &                                cl_parser%get('-siestaf'), flags,        &
     &                                cl_parser%get('-num_r', 1), num_p,       &
     &                                cl_parser%get('-num_z', 1),              &
     &                                cl_parser%get('-rmax', 0.0_dp),          &
     &                                cl_parser%get('-rmin', 0.0_dp),          &
     &                                cl_parser%get('-zmax', 0.0_dp),          &
     &                                cl_parser%get('-zmin', 0.0_dp),          &
     &                                parallel, io_unit)
      END IF
      CALL context%set_up_grid(cl_parser%get('-p_start', -1),                  &
     &                         cl_parser%get('-p_end', -1),                    &
     &                         parallel, io_unit)
      CALL context%write(cl_parser%get('-outf'), parallel)

      DEALLOCATE(context)
      DEALLOCATE(cl_parser)

      CALL profiler_set_stop_time('bmw_main', start_time)
      IF (parallel%offset .eq. 0) THEN
         CALL profiler_write(io_unit)
      END IF
      CALL profiler_destruct

      CLOSE(io_unit)

      DEALLOCATE(parallel)

1000  FORMAT('BMW ',i4,' Series.')
1001  FORMAT('Required flag ',a,' not set.')

      END PROGRAM

max_stars_repo_name: ORNL-Fusion/BMW

max_issues_repo_name: ORNL-Fusion/BMW

max_forks_repo_name: ORNL-Fusion/BMW


!F90 ISO_C_BINGING wrapper for socket communication.

!Copyright (C) 2013, Michele Ceriotti

!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.

!Contains both the functions that transmit data to the socket and read the data
!back out again once finished, and the function which opens the socket initially.

!Functions:
!   open_socket: Opens a socket with the required host server, socket type and
!      port number.
!   write_buffer: Writes a string to the socket.
!   read_buffer: Reads data from the socket.

   MODULE F90MESSAGES
   USE ISO_C_BINDING

   IMPLICIT NONE

  INTERFACE
     SUBROUTINE initialize_client_c() BIND(C, name="initialize_client")
       USE ISO_C_BINDING
     END SUBROUTINE initialize_client_c
  END INTERFACE

!    SUBROUTINE writebuffer_csocket(psockfd, pdata, plen) BIND(C, name="writebuffer")
!      USE ISO_C_BINDING
!    INTEGER(KIND=C_INT)                      :: psockfd
!    TYPE(C_PTR), VALUE                       :: pdata
!    INTEGER(KIND=C_INT)                      :: plen
!
!    END SUBROUTINE writebuffer_csocket

!    SUBROUTINE readbuffer_csocket(psockfd, pdata, plen) BIND(C, name="readbuffer")
!      USE ISO_C_BINDING
!    INTEGER(KIND=C_INT)                      :: psockfd
!    TYPE(C_PTR), VALUE                       :: pdata
!    INTEGER(KIND=C_INT)                      :: plen
!
!    END SUBROUTINE readbuffer_csocket
!  END INTERFACE

   CONTAINS

   SUBROUTINE initialize_client()
     CALL initialize_client_c()
   END SUBROUTINE

  END MODULE


text:
!******************************************************************************
! Copyright c 2013, The Ames Laboratory, Iowa State University, and Rutgers
! University*.  All rights reserved.
!
! This software was authored by Yongxin Yao, Nicola Lanata*, Gabriel Kotliar*,
! Cai-Zhuang Wang, and Kai-Ming Ho, at The Ames Laboratory and
! Rutgers University and was supported by the U.S.
! Department of Energy (DOE), Office of Science,
! Basic Energy Sciences, Materials Science and Engineering Division.
! The Ames Laboratory is operated by Iowa State University for DOE
! under U.S. Government contract DE-AC02-07CH11358.
! The U.S. Government has the rights to use, reproduce, and
! distribute this software.
! NEITHER THE GOVERNMENT, THE AMES LABORATORY, IOWA STATE UNIVERSITY,
! NOR RUTGERS UNIVERSITY MAKES ANY WARRANTY,
! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
! If software is modified to produce derivative works,
! such modified software should be clearly marked,
! so as to not confuse it with the version available from
! The Ames Laboratory and Rutgers University.
!
! Additionally, redistribution and use in source and binary forms,
! with or without modification,
! are permitted provided that the following conditions are met:
!
!     Redistribution of source code must retain the above copyright notice,
!     this list of conditions, and the following disclaimer.
!
!     Redistribution 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 distribution.
!
!     Neither the name of The Ames Laboratory, Iowa State University,
!     Rutgers University, the U.S. Government, nor the names of
!     its contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE AMES LABORATORY, IOWA STATE UNIVERSITY,
! RUTGERS UNIVERSITY, 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 GOVERNMENT, THE AMES LABORATORY,
! IOWA STATE UNIVERSITY, RUTGERS UNIVERSITY, OR 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.
!******************************************************************************

module sparse
    use gprec
    implicit none

    type dcsr_matrix
        integer :: nrow=0,ncol=0
        real(q),pointer :: a(:)=>null()
        integer,pointer :: i(:)=>null(),j(:)=>null()
    end type dcsr_matrix

    type zcsr_matrix
        integer :: nrow=0,ncol=0
        complex(q),pointer :: a(:)=>null()
        integer,pointer :: i(:)=>null(),j(:)=>null()
    end type zcsr_matrix

    type dcoo_matrix
        integer :: nrow=0,ncol=0,nnz=0
        real(q),pointer :: a(:)=>null()
        integer,pointer :: i(:)=>null(),j(:)=>null()
    end type dcoo_matrix

    type zvector
        complex(q),pointer :: a(:)=>null()
    end type zvector

    type ivector
        integer :: imin=0,imax=0
        integer,pointer :: i(:)=>null()
    end type ivector

    contains

    subroutine dcsr_syamuzx_sk(a,x,ax)
    type(dcsr_matrix),intent(in)::a
    complex(q),intent(in)::x(*)
    complex(q),intent(inout)::ax(*)

    integer i,k,j

    do i=1,a%nrow
        do k=a%i(i),a%i(i+1)-1
            j=a%j(k)
            ax(i)=ax(i)+a%a(k)*x(j)
            ! lower triangular part
            if(j==i)cycle
            ax(j)=ax(j)+a%a(k)*x(i)
        enddo
    enddo
    return

    end subroutine dcsr_syamuzx_sk

    subroutine dcsr_sylamuzx_sk(coef,a,x,ax)
    real(q),intent(in)::coef
    type(dcsr_matrix),intent(in)::a
    complex(q),intent(in)::x(*)
    complex(q),intent(inout)::ax(*)

    integer i,k,j

    do i=1,a%nrow
        do k=a%i(i),a%i(i+1)-1
            j=a%j(k)
            ax(i)=ax(i)+a%a(k)*x(j)*coef
            ! lower triangular part
            if(j==i)cycle
            ax(j)=ax(j)+a%a(k)*x(i)*coef
        enddo
    enddo
    return

    end subroutine dcsr_sylamuzx_sk

    subroutine zcsr_syamux_sk(a,x,ax)
    type(zcsr_matrix),intent(in)::a
    complex(q),intent(in)::x(*)
    complex(q),intent(inout)::ax(*)

    integer i,k,j

    do i=1,a%nrow
        do k=a%i(i),a%i(i+1)-1
            j=a%j(k)
            ax(i)=ax(i)+a%a(k)*x(j)
            ! lower triangular part
            if(j==i)cycle
            ax(j)=ax(j)+conjg(a%a(k))*x(i)
        enddo
    enddo
    return

    end subroutine zcsr_syamux_sk

    subroutine alloc_zcsr(a,nnz,nrow,ncol)
    integer nnz,nrow,ncol
    type(zcsr_matrix)::a

    integer ierr

    allocate(a%a(nnz),a%j(nnz),a%i(nrow+1),stat=ierr)
    if(ierr/=0)then
        write(0,'(" Error: allocation in alloc_zcsr with nnz = ",i0)')nnz
        stop 90
    endif
    a%nrow=nrow; a%ncol=ncol
    a%a=0; a%j=0; a%i=0
    return

    end subroutine alloc_zcsr

    subroutine copy_zcsr(a,b)
    type(zcsr_matrix)::a,b

    integer nnz

    nnz=a%i(a%nrow+1)-1
    b%nrow=a%nrow; b%ncol=a%ncol
    allocate(b%a(nnz),b%j(nnz),b%i(a%nrow+1))
    b%a=a%a(1:nnz)
    b%j=a%j(1:nnz)
    b%i=a%i
    return

    end subroutine copy_zcsr

    subroutine alloc_dcsr(a,nnz,nrow,ncol)
    integer nnz,nrow,ncol
    type(dcsr_matrix)::a

    allocate(a%a(nnz),a%j(nnz),a%i(nrow+1))
    a%nrow=nrow; a%ncol=ncol
    a%a=0; a%j=0; a%i=0
    return

    end subroutine alloc_dcsr

    subroutine alloc_dcoo(a,nnz,nrow,ncol)
    integer nnz,nrow,ncol
    type(dcoo_matrix)::a

    allocate(a%a(nnz),a%j(nnz),a%i(nnz))
    a%nrow=nrow; a%ncol=ncol; a%nnz=nnz
    a%a=0; a%j=0; a%i=0
    return

    end subroutine alloc_dcoo

    subroutine dealloc_zcsr(a)
    type(zcsr_matrix)::a

    deallocate(a%a,a%j,a%i)
    nullify(a%a,a%j,a%i)
    a%nrow=0; a%ncol=0
    return

    end subroutine dealloc_zcsr

    subroutine dealloc_dcoo(a)
    type(dcoo_matrix)::a

    deallocate(a%a,a%j,a%i)
    a%nrow=0; a%ncol=0; a%nnz=0
    return

    end subroutine dealloc_dcoo

    subroutine dealloc_dcsr(a)
    type(dcsr_matrix)::a

    deallocate(a%a,a%j,a%i)
    nullify(a%a,a%j,a%i)
    a%nrow=0; a%ncol=0
    return

    end subroutine dealloc_dcsr

    subroutine zs_dcoomux(zs,dcoo,x,y)
    type(dcoo_matrix),intent(in)::dcoo
    complex(q),intent(in)::zs,x(*)
    complex(q),intent(inout)::y(*)

    integer inz,i

    do inz=1,dcoo%nnz
        i=dcoo%i(inz)
        y(i)=y(i)+dcoo%a(inz)*x(dcoo%j(inz))*zs
    enddo
    return

    end subroutine zs_dcoomux

    subroutine zs_dcoohmux(zs,dcoo,x,y)
    type(dcoo_matrix),intent(in)::dcoo
    complex(q),intent(in)::zs,x(*)
    complex(q),intent(inout)::y(*)

    integer inz,j

    do inz=1,dcoo%nnz
        j=dcoo%j(inz)
        y(j)=y(j)+dcoo%a(inz)*x(dcoo%i(inz))*zs
    enddo
    return

    end subroutine zs_dcoohmux

    subroutine zvh_dcoo_zv(dcoo,x,y)
    type(dcoo_matrix),intent(in)::dcoo
    complex(q),intent(in)::x(*)
    complex(q),intent(inout)::y

    integer inz,i,j

    do inz=1,dcoo%nnz
        i=dcoo%i(inz)
        j=dcoo%j(inz)
        y=y+conjg(x(i))*dcoo%a(inz)*x(j)
    enddo
    return

    end subroutine zvh_dcoo_zv

    subroutine zv2h_sdcoo_zv1(zs,dcoo,x1,x2,y)
    type(dcoo_matrix),intent(in)::dcoo
    complex(q),intent(in)::zs,x1(*),x2(*)
    complex(q),intent(inout)::y

    integer inz,i,j

    do inz=1,dcoo%nnz
        i=dcoo%i(inz)
        j=dcoo%j(inz)
        y=y+zs*conjg(x2(i))*dcoo%a(inz)*x1(j)
    enddo
    return

    end subroutine zv2h_sdcoo_zv1

    subroutine zv2h_sdcooh_zv1(zs,dcoo,x1,x2,y)
    type(dcoo_matrix),intent(in)::dcoo
    complex(q),intent(in)::zs,x1(*),x2(*)
    complex(q),intent(inout)::y

    integer inz,i,j

    do inz=1,dcoo%nnz
        i=dcoo%i(inz)
        j=dcoo%j(inz)
        y=y+zs*conjg(x2(j))*dcoo%a(inz)*x1(i)
    enddo
    return

    end subroutine zv2h_sdcooh_zv1

    subroutine zvh_dcsr_zv(dcsr,x,y)
    type(dcsr_matrix),intent(in)::dcsr
    complex(q),intent(in)::x(*)
    complex(q),intent(inout)::y

    integer inz,i,j,k

    do i=1,dcsr%nrow; do k=dcsr%i(i),dcsr%i(i+1)-1
        j=dcsr%j(k)
        y=y+conjg(x(i))*dcsr%a(k)*x(j)
    enddo; enddo
    return

    end subroutine zvh_dcsr_zv

    subroutine zvh_sydcsr_zv(dcsr,x,y)
    type(dcsr_matrix),intent(in)::dcsr
    complex(q),intent(in)::x(*)
    complex(q),intent(inout)::y

    integer inz,i,j,k

    do i=1,dcsr%nrow; do k=dcsr%i(i),dcsr%i(i+1)-1
        j=dcsr%j(k)
        y=y+conjg(x(i))*dcsr%a(k)*x(j)
        if(i==j)cycle
        y=y+conjg(x(j))*dcsr%a(k)*x(i)
    enddo; enddo
    return

    end subroutine zvh_sydcsr_zv

    subroutine zvh_zcsr_zv(zcsr,x,y)
    type(zcsr_matrix),intent(in)::zcsr
    complex(q),intent(in)::x(*)
    complex(q),intent(inout)::y

    integer inz,i,j,k

    do i=1,zcsr%nrow; do k=zcsr%i(i),zcsr%i(i+1)-1
        j=zcsr%j(k)
        y=y+conjg(x(i))*zcsr%a(k)*x(j)
    enddo; enddo
    return

    end subroutine zvh_zcsr_zv

    ! zcsr = zs * dcsr
    subroutine zsdcsr(zs,dcsr,zcsr)
    complex(q),intent(in)::zs
    type(dcsr_matrix),intent(in)::dcsr
    type(zcsr_matrix),intent(out)::zcsr

    integer nnz

    nnz=dcsr%i(dcsr%nrow+1)-1
    call alloc_zcsr(zcsr,nnz,dcsr%nrow,dcsr%ncol)
    zcsr%a=zs*dcsr%a(1:nnz)
    zcsr%i=dcsr%i
    zcsr%j=dcsr%j(1:nnz)
    return

    end subroutine zsdcsr

    ! zcsrb = zs * zcsra
    subroutine zszcsr(zs,zcsra,zcsrb)
    complex(q),intent(in)::zs
    type(zcsr_matrix),intent(in)::zcsra
    type(zcsr_matrix),intent(out)::zcsrb

    integer nnz

    nnz=zcsra%i(zcsra%nrow+1)-1
    call alloc_zcsr(zcsrb,nnz,zcsra%nrow,zcsra%ncol)
    zcsrb%a=zs*zcsra%a(1:nnz)
    zcsrb%i=zcsra%i
    zcsrb%j=zcsra%j(1:nnz)
    return

    end subroutine zszcsr

    ! zsum = zcsr + zs * dcsr
    subroutine zcsrplsdcsr1(zs,dcsr,zcsr,zsum)
    complex(q),intent(in)::zs
    type(dcsr_matrix),intent(in)::dcsr
    type(zcsr_matrix),intent(in)::zcsr
    type(zcsr_matrix),intent(out)::zsum

    integer i,nnz,k,jpos
    integer iw(zcsr%ncol)

    nnz=dcsr%i(dcsr%nrow+1)+zcsr%i(zcsr%nrow+1)-2
    call alloc_zcsr(zsum,nnz,zcsr%nrow,zcsr%ncol)
    zsum%i(1)=1
    nnz=0
    do i=1,zcsr%nrow
        iw=0
        do k=zcsr%i(i),zcsr%i(i+1)-1
            nnz=nnz+1
            zsum%j(nnz)=zcsr%j(k)
            zsum%a(nnz)=zcsr%a(k)
            iw(zcsr%j(k))=nnz
        enddo
        do k=dcsr%i(i),dcsr%i(i+1)-1
            jpos=iw(dcsr%j(k))
            if(jpos==0)then
                nnz=nnz+1
                zsum%j(nnz)=dcsr%j(k)
                zsum%a(nnz)=zs*dcsr%a(k)
                iw(dcsr%j(k))=nnz
            else
                zsum%a(jpos)=zsum%a(jpos)+zs*dcsr%a(k)
            endif
        enddo
        zsum%i(i+1)=nnz+1
    enddo
    return

    end subroutine zcsrplsdcsr1

    ! zcsr += zs * dcsr
    subroutine zcsrplsdcsr(zs,dcsr,zcsr)
    complex(q),intent(in)::zs
    type(dcsr_matrix),intent(in)::dcsr
    type(zcsr_matrix),intent(inout)::zcsr

    type(zcsr_matrix)::zsum

    if(zcsr%nrow==0)then
        call zsdcsr(zs,dcsr,zcsr)
    else
        call zcsrplsdcsr1(zs,dcsr,zcsr,zsum)
        call dealloc_zcsr(zcsr)
        call copy_zcsr(zsum,zcsr)
        call dealloc_zcsr(zsum)
    endif
    return

    end subroutine zcsrplsdcsr

    ! zsum = zs*zcsra + zcsrb
    subroutine zcsraplszcsrb1(zs,zcsra,zcsrb,zsum)
    complex(q),intent(in)::zs
    type(zcsr_matrix),intent(in)::zcsra,zcsrb
    type(zcsr_matrix),intent(out)::zsum

    integer i,nnz,k,jpos
    integer iw(zcsra%ncol)

    nnz=zcsra%i(zcsra%nrow+1)+zcsrb%i(zcsrb%nrow+1)-2
    call alloc_zcsr(zsum,nnz,zcsra%nrow,zcsra%ncol)
    zsum%i(1)=1
    nnz=0
    do i=1,zcsrb%nrow
        iw=0
        do k=zcsrb%i(i),zcsrb%i(i+1)-1
            nnz=nnz+1
            zsum%j(nnz)=zcsrb%j(k)
            zsum%a(nnz)=zcsrb%a(k)
            iw(zcsrb%j(k))=nnz
        enddo
        do k=zcsra%i(i),zcsra%i(i+1)-1
            jpos=iw(zcsra%j(k))
            if(jpos==0)then
                nnz=nnz+1
                zsum%j(nnz)=zcsra%j(k)
                zsum%a(nnz)=zs*zcsra%a(k)
                iw(zcsra%j(k))=nnz
            else
                zsum%a(jpos)=zsum%a(jpos)+zs*zcsra%a(k)
            endif
        enddo
        zsum%i(i+1)=nnz+1
    enddo
    return

    end subroutine zcsraplszcsrb1

    ! zsum += zs * zcsr
    subroutine zcsrplszcsra(zs,zcsr,zsum)
    complex(q),intent(in)::zs
    type(zcsr_matrix),intent(in)::zcsr
    type(zcsr_matrix),intent(inout)::zsum

    type(zcsr_matrix)::zsum_

    if(zsum%nrow==0)then
        call zszcsr(zs,zcsr,zsum)
    else
        call zcsraplszcsrb1(zs,zcsr,zsum,zsum_)
        call dealloc_zcsr(zsum)
        call copy_zcsr(zsum_,zsum)
        call dealloc_zcsr(zsum_)
    endif
    return

    end subroutine zcsrplszcsra

    subroutine calc_nnz_csrab(nrow,ncol,ia,ja,ib,jb,nnz)
    integer,intent(in)::nrow,ncol,ia(*),ja(*),ib(*),jb(*)
    integer,intent(out)::nnz

    integer i,j,j1,k,k1
    logical lzero(ncol)

    nnz=0
    do i=1,nrow
        lzero=.true.
        do j1=ia(i),ia(i+1)-1
            j=ja(j1)
            do k1=ib(j),ib(j+1)-1
                k=jb(k1)
                if(lzero(k))then
                    nnz=nnz+1
                    lzero(k)=.false.
                endif
            enddo
        enddo
    enddo
    return

    end subroutine calc_nnz_csrab

    ! zcsrab = zcsra * zcsrb
    subroutine zcsramuzcsrb1(zcsra,zcsrb,zcsrab)
    type(zcsr_matrix),intent(in)::zcsra,zcsrb
    type(zcsr_matrix),intent(out)::zcsrab

    integer i,nnz,j1,j,k1,k,jpos
    integer iw(zcsrb%ncol)

    call calc_nnz_csrab(zcsra%nrow,zcsrb%ncol,zcsra%i,zcsra%j,zcsrb%i, &
            &zcsrb%j,nnz)
    call alloc_zcsr(zcsrab,nnz,zcsra%nrow,zcsrb%ncol)
    zcsrab%i(1)=1
    nnz=0
    do i=1,zcsra%nrow
        iw=0
        do j1=zcsra%i(i),zcsra%i(i+1)-1
            j=zcsra%j(j1)
            do k1=zcsrb%i(j),zcsrb%i(j+1)-1
                k=zcsrb%j(k1)
                jpos=iw(k)
                if(jpos==0)then
                    nnz=nnz+1
                    zcsrab%j(nnz)=k
                    zcsrab%a(nnz)=zcsra%a(j1)*zcsrb%a(k1)
                    iw(k)=nnz
                else
                    zcsrab%a(jpos)=zcsrab%a(jpos)+zcsra%a(j1)*zcsrb%a(k1)
                endif
            enddo
        enddo
        zcsrab%i(i+1)=nnz+1
    enddo
    return

    end subroutine zcsramuzcsrb1

    ! zsum += s * zcsra * zcsrb
    subroutine zcsradd_szcsramuzcsrb(zs,zcsra,zcsrb,zsum)
    complex(q),intent(in)::zs
    type(zcsr_matrix),intent(in)::zcsra,zcsrb
    type(zcsr_matrix),intent(inout)::zsum

    type(zcsr_matrix)::zcsrab

    if(zsum%nrow==0)then
        call zcsramuzcsrb1(zcsra,zcsrb,zsum)
        if(abs(zs-1._q)>1.d-12)then
            zsum%a=zs*zsum%a
        endif
    else
        call zcsramuzcsrb1(zcsra,zcsrb,zcsrab)
        call zcsrplszcsra(zs,zcsrab,zsum)
        call dealloc_zcsr(zcsrab)
    endif
    return

    end subroutine zcsradd_szcsramuzcsrb

    ! zcsr = s_{a,b}*(dcsra(a,b) + dcsrb(a,b))
    subroutine zsdcsr_22darray_sum(n,s,dcsra,dcsrb,zcsr)
    integer,intent(in)::n
    complex(q),intent(in)::s(n,n)
    type(dcsr_matrix),intent(in)::dcsra(n,n),dcsrb(n,n)
    type(zcsr_matrix),intent(out)::zcsr

    integer i,j

    do i=1,n; do j=1,n
        if(abs(s(i,j))<1.d-12)cycle
        call zcsrplsdcsr(s(i,j),dcsra(i,j),zcsr)
        call zcsrplsdcsr(s(i,j),dcsrb(i,j),zcsr)
    enddo; enddo
    return

    end subroutine zsdcsr_22darray_sum

    subroutine zcsr_array_sqsum(n,zcsr,zsum)
    integer,intent(in)::n
    type(zcsr_matrix),intent(in)::zcsr(n)
    type(zcsr_matrix),intent(out)::zsum

    integer i
    complex(q),parameter::z1=cmplx(1._q,0._q)

    do i=1,n
        call zcsradd_szcsramuzcsrb(z1,zcsr(i),zcsr(i),zsum)
    enddo
    return

    end subroutine zcsr_array_sqsum

end module sparse

max_stars_repo_name: comscope/comsuite

max_issues_repo_name: comscope/Comsuite

max_forks_repo_name: comscope/Comsuite


text:
      subroutine dimdro(d,j)
C.....subroutine dimdro  01/06/75
      common/gfgel/gf(28)
      common/gftab/irdb(5),rdb(17)
      common/gfbftr/xtr,ytr
      dimension buf(28)
      f=rdb(17)/rdb(1)
      call wrstr(buf)
      m1=1
      m2=1
      if(gf(9).ge.0.) call strmod(0)
      i=abs(gf(9))
      if(iabs(j).eq.10) goto 15
      if(j) 16,15,17
 16   m1=0
      go to 15
 17   m2=0
 15   l=gf(i)
      if(gf(9)) 41,3,38
  38  call strmod(0)
  41  if(l) 10,4,8
   4  call rdmva(x0,y0,dl,th)
      call move(x0,y0,0)
      if(d.eq.0.)go to 6
      if(th-90.) 11,11,12
  11  th1=th+sign(90.,d)
      m3=1
      goto 7
  12  if(th.gt.270.) goto 11
      th1=th-sign(90.,d)
      m3=0
  7    call mova(abs(d)+0.3*f,th1,m1*m3+m2*(1-m3))
       call mova(dl,th,0)
       call mova(abs(d)+0.3*f,th1+180.,m1*(1-m3)+m2*m3)
       call mova(abs(d),th1,0)
   20  if(iabs(j).ge.10) goto 19
       call wrmva(dl,th+180.)
       call arrow(0)
   40  call move(x0,y0,0)
       call mova(dl,th,0)
      goto 9
   6  call mova(dl,th,0)
      goto 20
  19  call mova(0.6*f,th,0)
      call wrmva(0.6*f,th+180.)
      call arrow(1)
      call mova(dl,th+180.,1)
      call wrmva(0.6*f,th+180.)
      call arrow(-1)
      goto 40
  8   call rdaca(x0,y0,r,th0,thf)
      call rdacb(x0,y0,r,xf,yf,j1)
      call move(x0,y0,0)
      if(d.eq.0.)goto 22
      if(r*(thf-th0)) 23,3,24
  23  th1=180.+th0/2+thf/2
      goto 25
  24  th1=th0/2+thf/2
 25   k=0
      if(th1.gt.180.) k=1
  28  if(r*k) 29,30,31
  30  if(r) 31,3,29
  29  m3=0
      goto 32
  31  m3=1
  32  if(d*(0.5-k)) 34,3,33
  33  q=0.
      goto 35
  34  q=1.
  35  call mova(abs(d)+0.3*f,th0+q*180.,m1*m3+m2*(1-m3))
      call move(xf,yf,0)
      call mova(abs(d)+0.3*f,thf+q*180.,m1*(1-m3)+m2*m3)
      call mova(0.3*f,thf+q*180.+180.,0)
  39  xf1=xtr
      yf1=ytr
      r1=abs(r)+abs(d)*(1.-2.*q)
      if(iabs(j).ge.10) goto 37
      call wraca(sign(r1,-r),thf,th0)
      call arrow(0)
      call move(xf,yf,0)
      goto 9
  22  call move(xf,yf,0)
      goto 39
  37  call wraca(sign(r1,r),thf,thf+sign(8.,r))
      call arrow(-1)
      call move(xf1,yf1,0)
      call arcia(sign(r1,-r),thf,th0)
      call wraca(sign(r1,-r),th0,th0-sign(8.,r))
      call arrow(-1)
      call move(xf,yf,0)
      goto 9
  10  call gffals(6)
   9  call rdstr(buf)
   3  return
      end

ext: f

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
*> \brief \b DCHKLQT
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DCHKLQT( THRESH, TSTERR, NM, MVAL, NN, NVAL, NNB,
*                           NBVAL, NOUT )
*
*       .. Scalar Arguments ..
*       LOGICAL            TSTERR
*       INTEGER            NM, NN, NNB, NOUT
*       DOUBLE PRECISION   THRESH
*       ..
*       .. Array Arguments ..
*       INTEGER            MVAL( * ), NBVAL( * ), NVAL( * )
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DCHKLQT tests DGELQT and DGEMLQT.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] THRESH
*> \verbatim
*>          THRESH is DOUBLE PRECISION
*>          The threshold value for the test ratios.  A result is
*>          included in the output file if RESULT >= THRESH.  To have
*>          every test ratio printed, use THRESH = 0.
*> \endverbatim
*>
*> \param[in] TSTERR
*> \verbatim
*>          TSTERR is LOGICAL
*>          Flag that indicates whether error exits are to be tested.
*> \endverbatim
*>
*> \param[in] NM
*> \verbatim
*>          NM is INTEGER
*>          The number of values of M contained in the vector MVAL.
*> \endverbatim
*>
*> \param[in] MVAL
*> \verbatim
*>          MVAL is INTEGER array, dimension (NM)
*>          The values of the matrix row dimension M.
*> \endverbatim
*>
*> \param[in] NN
*> \verbatim
*>          NN is INTEGER
*>          The number of values of N contained in the vector NVAL.
*> \endverbatim
*>
*> \param[in] NVAL
*> \verbatim
*>          NVAL is INTEGER array, dimension (NN)
*>          The values of the matrix column dimension N.
*> \endverbatim
*>
*> \param[in] NNB
*> \verbatim
*>          NNB is INTEGER
*>          The number of values of NB contained in the vector NBVAL.
*> \endverbatim
*>
*> \param[in] NBVAL
*> \verbatim
*>          NBVAL is INTEGER array, dimension (NBVAL)
*>          The values of the blocksize NB.
*> \endverbatim
*>
*> \param[in] NOUT
*> \verbatim
*>          NOUT is INTEGER
*>          The unit number for output.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_lin
*
*  =====================================================================
      SUBROUTINE DCHKLQT( THRESH, TSTERR, NM, MVAL, NN, NVAL, NNB,
     $                     NBVAL, NOUT )
      IMPLICIT NONE
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            TSTERR
      INTEGER            NM, NN, NNB, NOUT
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      INTEGER            MVAL( * ), NBVAL( * ), NVAL( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            NTESTS
      PARAMETER          ( NTESTS = 6 )
*     ..
*     .. Local Scalars ..
      CHARACTER*3        PATH
      INTEGER            I, J, K, T, M, N, NB, NFAIL, NERRS, NRUN,
     $                   MINMN
*
*     .. Local Arrays ..
      DOUBLE PRECISION   RESULT( NTESTS )
*     ..
*     .. External Subroutines ..
      EXTERNAL           ALAERH, ALAHD, ALASUM, DERRLQT, DLQT04
*     ..
*     .. Scalars in Common ..
      LOGICAL            LERR, OK
      CHARACTER*32       SRNAMT
      INTEGER            INFOT, NUNIT
*     ..
*     .. Common blocks ..
      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
      COMMON             / SRNAMC / SRNAMT
*     ..
*     .. Executable Statements ..
*
*     Initialize constants
*
      PATH( 1: 1 ) = 'D'
      PATH( 2: 3 ) = 'TQ'
      NRUN = 0
      NFAIL = 0
      NERRS = 0
*
*     Test the error exits
*
      IF( TSTERR ) CALL DERRLQT( PATH, NOUT )
      INFOT = 0
*
*     Do for each value of M in MVAL.
*
      DO I = 1, NM
         M = MVAL( I )
*
*        Do for each value of N in NVAL.
*
         DO J = 1, NN
            N = NVAL( J )
*
*        Do for each possible value of NB
*
            MINMN = MIN( M, N )
            DO K = 1, NNB
               NB = NBVAL( K )
*
*              Test DGELQT and DGEMLQT
*
               IF( (NB.LE.MINMN).AND.(NB.GT.0) ) THEN
                  CALL DLQT04( M, N, NB, RESULT )
*
*                 Print information about the tests that did not
*                 pass the threshold.
*
                  DO T = 1, NTESTS
                     IF( RESULT( T ).GE.THRESH ) THEN
                        IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                       CALL ALAHD( NOUT, PATH )
                        WRITE( NOUT, FMT = 9999 )M, N, NB,
     $                       T, RESULT( T )
                        NFAIL = NFAIL + 1
                     END IF
                  END DO
                  NRUN = NRUN + NTESTS
               END IF
            END DO
         END DO
      END DO
*
*     Print a summary of the results.
*
      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
*
 9999 FORMAT( ' M=', I5, ', N=', I5, ', NB=', I4,
     $      ' test(', I2, ')=', G12.5 )
      RETURN
*
*     End of DCHKLQT
*
      END

max_stars_repo_name: sergey-v-kuznetsov/lapack

max_issues_repo_name: sergey-v-kuznetsov/lapack

max_forks_repo_name: sergey-v-kuznetsov/lapack


text:
Amy Googenspa once bought a pint of Ben & Jerrys and ate the entire thing.

20070102 20:42:37 nbsp Sorry, but no.  Not even really a fan of their food.  I just call people on idiotic comments like Crepeville only hires Mexicans or I had to drink out of a plastic cup, boo hoo! Users/JeffSpeckles

20070806 19:50:39 nbsp I noticed your comment on the Blockbuster page and wondered if you ever filed a complaint.  If you filed it with the management (esp. in the Downtown store) know that it didnt get there because when I worked there she couldnt have cared less what customers thought.  I urge you (and anyone else you know who has complaints) to call the corporation.  I doubt theyll do anything (they treat their employees as bad as their customers) but its worth a shot.  I tried really hard when I was employed to treat customers right but was actively discouraged because it often took too long.  Thanks, from a former (and very disgruntled) employee Users/Naomi

20080317 18:45:24 nbsp {{{Agricultural and Resource Economics/Talk}}} Users/JabberWokky

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


REAL FUNCTION addfloat( A, B )
      IMPLICIT NONE
      REAL :: A, B

      addfloat = A + B

END FUNCTION addfloat

DOUBLE PRECISION FUNCTION adddouble( A, B )
      IMPLICIT NONE
      DOUBLE PRECISION :: A, B

      adddouble = A + B

END FUNCTION adddouble