text:
      SUBROUTINE PSTQ1(NTYPE)
C  THIS ROUTINE CALCULATES PHASE I OUTPUT FOR PLA3
C  FOR COMBINATION ELEMENTS
C
C**************** PHASE I  STRESS DATA RECOVERY ************************
C **********************************************************************
C
C     9/12/67         E C P T     L I S T I N G
C                    ***************************
C ECPT  TRMEM   QDMEM   TRPLT   QDPLT   TRIA1   QUAD1   TRIA2   QUAD2
C **********************************************************************
C   1   EL.ID   EL.ID   EL.ID   EL.ID   EL.ID   EL.ID   EL.ID   EL.ID
C   2   GRID A  GRID A  GRID A  GRID A  GRID A  GRID A  GRID A  GRID A
C   3   GRID B  GRID B  GRID B  GRID B  GRID B  GRID B  GRID B  GRID B
C   4   GRID C  GRID C  GRID C  GRID C  GRID C  GRID C  GRID C  GRID C
C   5   THETA   GRID D  THETA   GRID D  THETA   GRID D  THETA   GRID D
C   6   MATID   THETA   MATID1  THETA   MATID1  THETA   MAT ID  THETA
C   7   T       MAT ID  I       MATID1  T1      MATID1  T       MAT ID
C   8   NS MASS T       MATID2  I       MATID2  T1      NS MASS T
C   9   CSID 1  NS MASS T2      MATID2  I       MATID2  CSID 1  NS MASS
C  10   X1      CSID 1  NS MASS T2      MATID3  I       X1      CSID 1
C  11   Y1      X1      Z1      NS MASS T2      MATID3  Y1      X1
C  12   Z1      Y1      Z2      Z1      NS MASS T2      Z1      Y1
C  13   CSID 2  Z1      CSID 1  Z2      Z1      NS MASS CSID 2  Z1
C  14   X2      CSID 2  X1      CSID 1  Z2      Z1      X2      CSID 2
C  15   Y2      X2      Y1      X1      CSID 1  Z2      Y2      X2
C  16   Z2      Y2      Z1      Y1      X1      CSID 1  Z2      Y2
C  17   CSID 3  Z2      CSID 2  Z1      Y1      X1      CSID 3  Z2
C  18   X3      CSID 3  X2      CSID 2  Z1      Y1      X3      CSID 3
C  19   Y3      X3      Y2      X2      CSID 2  Z1      Y3      X3
C  20   Z3      Y3      Z2      Y2      X2      CSID 2  Z3      Y3
C  21   TEMP    Z3      CSID 3  Z2      Y2      X2      TEMP    Z3
C  22           CSID 4  X3      CSID 3  Z2      Y2              CSID 4
C  23           X4      Y3      X3      CSID 3  Z2              X4
C  24           Y4      Z3      Y3      X3      CSID 3          Y4
C  25           Z4      TEMP    Z3      Y3      X3              Z4
C  26           TEMP            CSID 4  Z3      Y3              TEMP
C  27                           X4      TEMP    Z3
C  28                           Y4              CSID 4
C  29                           Z4              X4
C  30                           TEMP            Y4
C  31                                           Z4
C  32                                           TEMP
C **********************************************************************
C
      DIMENSION SAVE(32)
C
      COMMON /PLA3ES/ ECPT(100), PH1OUT(173) ,DUMMY(27)
C
C
C     THIS SUBROUTINE INCORPORATES TRIA1, QUAD1, TRIA2, QUAD2
C
C              NTYPE = 1  IMPLIES STRIA1
C              NTYPE = 2  IMPLIES STRIA2
C              NTYPE = 3  IMPLIES SQUAD1
C              NTYPE = 4  IMPLIES SQUAD2
C
C     SAVE THE INCOMING ECPT
C
      DO 10 I=1,32
   10 SAVE(I) = ECPT(I)
C
C     TRANSFER TO OPERATIONS DESIRED
C
C              STRIA1    STRIA2    SQUAD1    SQUAD2
      GO TO(20,100,150,230),NTYPE
C
C     **************
C     *** STRIA1 ***
C     **************
C
C     SET UP ECPT FOR PSTRM1, FIRST CHECK T1 FOR ZERO
   20 IF( SAVE(7) .EQ. 0.0E0 ) GO TO 50
      DO 30 I=9,21
   30 ECPT(I) = SAVE(I + 6)
C
      CALL PSTRM1 (0)
C
C     MOVE OUTPUT FROM PSTRM1 TO NEAR BOTTOM OF PH1OUT
C     WORDS (1 THRU 36) DOWN TO (99 THRU 134)
C
C
      DO 40 I=1,36
   40 PH1OUT(I + 98) = PH1OUT(I)
      GO TO 60
C
   50 PH1OUT( 99) = ECPT(1)
      PH1OUT(100) = 0.0E0
C
C     SET UP CALL TO PSTPL1, CHECK I EQUAL TO ZERO
   60 IF( SAVE(9) .EQ. 0.0E0 ) GO TO 90
      DO 70 I=1,5
   70 ECPT(I) = SAVE(I)
      DO 80 I=6,25
   80 ECPT(I) = SAVE(I + 2)
C
      CALL PSTPL1
      RETURN
C
   90 PH1OUT(1) = ECPT(1)
      PH1OUT(2) = 0.0E0
      RETURN
C
C     **************
C     *** STRIA2 ***
C     **************
  100 IF( SAVE(7) .EQ. 0.0E0 ) GO TO 140
C     SET UP CALL TO PSTRM1
C
C      ECPT IS OK AS DELIVERED TO THIS ROUTINE
C
      CALL PSTRM1(0)
C
C     MOVE OUTPUT FROM PSTRM1 TO NEAR BOTTOM OF PH1OUT
C     WORDS (1 THRU 36) DOWN TO (99 THRU 134)
C
      DO 110 I=1,36
  110 PH1OUT(I + 98) = PH1OUT(I)
C
C     SET UP CALL TO PSTPL1
C
      DO 120 I=1,6
  120 ECPT(I) = SAVE(I)
      ECPT(7) = SAVE(7) ** 3  / 12.0E0
      ECPT(8) = SAVE(6)
      ECPT(9) = SAVE(7)
      ECPT(10)= SAVE(8)
      ECPT(11) = SAVE(7)/2.0E0
      ECPT(12) = -ECPT(11)
      DO 130 I=13,25
  130 ECPT(I) = SAVE(I - 4)
C
      CALL PSTPL1
      RETURN
C
  140 PH1OUT(  1) = ECPT(1)
      PH1OUT(  2) = 0.0E0
      PH1OUT( 99) = ECPT(1)
      PH1OUT(100) = 0.0E0
      RETURN
C
C     **************
C     *** SQUAD1 ***
C     **************
C
  150 IF(SAVE(8).EQ.0.0E0)GO TO 180
C
C     SET UP CALL TO PSQDM1
C
      ECPT(9) = SAVE(13)
      DO 160 I=10,26
  160 ECPT(I) = SAVE(I+6)
C
      CALL PSQDM1
C
C     MOVE OUTPUT DOWN TO NEAR BOTTOM OF PH1OUT
C     WORDS (1 THRU 45) DOWN TO (129 THRU 173)
C
      DO 170 I=1,45
  170 PH1OUT(I + 128) = PH1OUT(I)
C
      GO TO 190
  180 PH1OUT(129) = ECPT(1)
      PH1OUT(130) = 0.0E0
C
  190 IF( SAVE(10) .EQ. 0.0E0 ) GO TO 220
C
C     SET UP CALL TO PSQPL1
C
      DO 200 I=1,6
  200 ECPT(I) = SAVE(I)
      DO 210 I=7,30
  210 ECPT(I) = SAVE(I + 2)
C
      CALL PSQPL1
      RETURN
C
  220 PH1OUT(1) = ECPT(1)
      PH1OUT(2) = 0.0E0
      RETURN
C
C     **************
C     *** SQUAD2 ***
C     **************
C
  230 IF( SAVE(8) .EQ. 0.0E0 ) GO TO 270
C
C     SET UP CALL TO PSQDM1
C
C      ECPT IS OK AS DELIVERED TO THIS ROUTINE
C
      CALL PSQDM1
C
C     MOVE OUTPUT DOWN TO NEAR BOTTOM OF PH1OUT
C     WORDS (1 THRU 45) DOWN TO (129 THRU 173)
C
      DO 240 I=1,45
  240 PH1OUT(I + 128) = PH1OUT(I)
C
C
C     SET UP CALL TO PSQPL1
C
      DO 250 I=1,7
  250 ECPT(I) = SAVE(I)
      ECPT(8) = SAVE(8) **3 / 12.0E0
      ECPT(9) = SAVE(7)
      ECPT(10)= SAVE(8)
      ECPT(11)= SAVE(9)
      ECPT(12) = SAVE(8)/2.0E0
      ECPT(13) =-ECPT(12)
      DO 260 I=14,30
  260 ECPT(I) = SAVE(I - 4)
C
      CALL PSQPL1
C
      RETURN
C
  270 PH1OUT(1) = ECPT(1)
      PH1OUT(2) = 0.0E0
      PH1OUT(129) = ECPT(1)
      PH1OUT(130) = 0.0E0
      RETURN
      END


text:
*----------------------------------------------------------------------*
      subroutine diag_guess(mel_out,xlist,idxlist,nlist,iroot,isign,
     &     op_info,str_info,strmap_info,orb_info)
*----------------------------------------------------------------------*
*     initialize list with diagonal guess for (iroot)th root
*     using the arrays xlist, idxlist which provide the lowest entries
*     (+ indices) of some diagonal approximation
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_orbinf.h'
      include 'def_graph.h'
      include 'mdef_operator_info.h'
      include 'def_strinf.h'
      include 'def_strmapinf.h'
      include 'ifc_memman.h'

      real(8), parameter ::
     &     thresh = 1d-10

      type(me_list), intent(inout) ::
     &     mel_out
      integer, intent(in) ::
     &     nlist, isign, iroot
      integer, intent(in) ::
     &     idxlist(nlist)
      real(8), intent(in) ::
     &     xlist(nlist)
      type(orbinf), intent(in) ::
     &     orb_info
      type(strinf), intent(in) ::
     &     str_info
      type(strmapinf), intent(in) ::
     &     strmap_info
      type(operator_info), intent(in) ::
     &     op_info

      integer ::
     &     jroot, idx, nset
      integer ::
     &     idxset(2)
      real(8) ::
     &     valset(2)

      integer ::
     &     idxlist_ba(nlist)

      ! if requested, get alpha<->beta interchanged indices
      if (isign.ne.0)
     &     call symidx_ab(idxlist_ba,
     &     idxlist,nlist,mel_out,
     &     op_info,str_info,strmap_info,orb_info)

      ! a) find out which elements to set
      jroot = 0
      idx = 1
      nset = 0
      do idx = 1, nlist
        if (isign.ne.0) then
          if (idxlist(idx).lt.abs(idxlist_ba(idx))) then
            ! set a spin combination
            jroot = jroot+1
            nset = 2
            idxset(1) = abs(idxlist(idx))
            idxset(2) = abs(idxlist_ba(idx))
            valset(1) = 1d0/sqrt(2d0)
            valset(2) = dble(isign)*dble(sign(1,idxlist_ba(idx)))
     &                  /sqrt(2d0)
          else if (idxlist(idx).eq.abs(idxlist_ba(idx))
     &           .and.isign.eq.+1) then
            if (idxlist_ba(idx).lt.0)
     &           call quit(1,'diag_guess','unexpected case')
            ! set a single element
            jroot = jroot+1
            nset = 1
            idxset(1) = abs(idxlist(idx))
            valset(1) = 1d0
          end if
        else
          ! set a single element
          jroot = jroot+1
          nset = 1
          idxset(1) = abs(idxlist(idx))
          valset(1) = 1d0
        end if
        if (jroot.eq.iroot) exit
      end do

      if (nset.eq.0)
     &     call quit(1,'diag_guess','nset==0 ??')

      if (iprlvl.ge.5) then
        if (nset.eq.1) then
          write(lulog,'(3x,"guess for root ",i3,": ",f12.6,i10)')
     &         iroot,valset(1),idxset(1)
        else
          write(lulog,'(3x,"guess for root ",i3,": ",2(f12.6,i10))')
     &         iroot,valset(1),idxset(1),valset(2),idxset(2)
        end if
      end if

      call set_list(mel_out,idxset,valset,nset)

      return
      end

ext: f

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
Christopher Tuffley was a master mathematician and part of the Math math department at UC Davis, through Fall 2006. He taught Math 22A, Linear Algebra, and Math 25, advanced calculus, during Fall of 2006, and was known for his great teaching style. A few random facts is that he is from Australia (hence zed not zee and diagonalisation), his aides never showed up except for tests, he had a mysterious bump on his forehead, and lastly he can juggle quite well as he demonstrated before his Fall 2006 Math 22A final examination.

As his UCD math department http://math.ucdavis.edu/~tuffley website indicates, Dr. Tuffley is now teaching at Massey University, in New Zealand.

   Would it simplify the Mathematics page to just put all of the above on that site, under Christopher Tuffley? That, rather than have two instructors w/ bulleted info, and one w/ a linked subpage (this one)? Users/LeonardMarque
     Probably, but Id imagine the eventual goal is to give each professor their own entry.  Id say being a professor makes you a public figure, certainly for the student segment of wiki users (and usually for everybody, since they have usually published a book or some research of note). Users/JabberWokky jw
     I agree with Leonard. Since this guy moved away, there is probably not going to be much more to say about him.  Users/NickSchmalenberger

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------

module m_hecmw_matrix_ordering_MC
  use hecmw_util
  implicit none

  private
  public :: hecmw_matrix_ordering_MC

contains

  subroutine hecmw_matrix_ordering_MC(N, indexL, itemL, indexU, itemU, &
      perm_cur, ncolor_in, ncolor_out, COLORindex, perm, iperm)
    implicit none
    integer(kind=kint), intent(in) :: N
    integer(kind=kint), intent(in) :: indexL(0:), indexU(0:)
    integer(kind=kint), intent(in) :: itemL(:), itemU(:)
    integer(kind=kint), intent(in) :: perm_cur(:)
    integer(kind=kint), intent(in) :: ncolor_in
    integer(kind=kint), intent(out) :: ncolor_out
    integer(kind=kint), intent(out) :: COLORindex(0:)
    integer(kind=kint), intent(out) :: perm(:), iperm(:)
    integer(kind=kint), allocatable :: iwk(:)
    integer(kind=kint) :: nn_color, cntall, cnt, color
    integer(kind=kint) :: i, inode, j, jnode
    allocate(iwk(N))
    iwk = 0
    nn_color = N / ncolor_in
    cntall = 0
    COLORindex(0) = 0
    do color=1,N
      cnt = 0
      do i=1,N
        inode = perm_cur(i)
        if (iwk(inode) > 0 .or. iwk(inode) == -1) cycle
        ! if (iwk(inode) == 0)
        iwk(inode) = color
        cntall = cntall + 1
        perm(cntall) = inode
        cnt = cnt + 1
        if (cnt == nn_color) exit
        if (cntall == N) exit
        ! mark all connected and uncolored nodes
        do j = indexL(inode-1)+1, indexL(inode)
          jnode = itemL(j)
          if (iwk(jnode) == 0) iwk(jnode) = -1
        end do
        do j = indexU(inode-1)+1, indexU(inode)
          jnode = itemU(j)
          if (jnode > N) cycle
          if (iwk(jnode) == 0) iwk(jnode) = -1
        end do
      end do
      COLORindex(color) = cntall
      if (cntall == N) then
        ncolor_out = color
        exit
      end if
      ! unmark all marked nodes
      do i=1,N
        if (iwk(i) == -1) iwk(i) = 0
      end do
    end do
    deallocate(iwk)
    ! make iperm
    do i=1,N
      iperm(perm(i)) = i
    end do
  end subroutine hecmw_matrix_ordering_MC

end module m_hecmw_matrix_ordering_MC


text:
../../src/vlog/async_fifo.v
../../src/vlog/fifo_2mem.v
../../src/vlog/rptr_empty.v
../../src/vlog/sync_r2w.v
../../src/vlog/sync_w2r.v
../../src/vlog/wptr_full.v


text:
C
C	$Id: cshval.f,v 1.5 2008-07-27 03:10:07 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      DOUBLE PRECISION FUNCTION CSHVAL (B,H1,H2,HP1,HP2,SIGMA)
      DOUBLE PRECISION B, H1, H2, HP1, HP2, SIGMA
C
C***********************************************************
C
C                                              From SSRFPACK
C                                            Robert J. Renka
C                                  Dept. of Computer Science
C                                       Univ. of North Texas
C                                           renka@cs.unt.edu
C                                                   11/21/96
C
C   Given a line segment P1-P2 containing a point P, along
C with values and derivatives at the endpoints, this func-
C tion returns the value H(P), where H is the Hermite inter-
C polatory tension spline defined by the endpoint data.
C
C On input:
C
C       B = Local coordinate of P with respect to P1-P2:
C           P = B*P1 + (1-B)*P2, and thus B = d(P,P2)/
C           d(P1,P2), where d(P1,P2) is the distance between
C           P1 and P2.  B < 0 or B > 1 results in extrapola-
C           tion.
C
C       H1,H2 = Values interpolated at P1 and P2, respec-
C               tively.
C
C       HP1,HP2 = Products of d(P1,P2) with first order der-
C                 ivatives at P1 and P2, respectively.  HP1
C                 may, for example, be the scalar product of
C                 P2-P1 with a gradient at P1.
C
C       SIGMA = Nonnegative tension factor associated with
C               the spline.  SIGMA = 0 corresponds to a
C               cubic spline, and H approaches the linear
C               interpolant of H1 and H2 as SIGMA increases.
C
C Input parameters are not altered by this function.
C
C On output:
C
C       CSHVAL = Interpolated value H(P).
C
C SSRFPACK module required by CSHVAL:  CSSNHCSH
C
C Intrinsic functions called by CSHVAL:  ABS, EXP
C
C***********************************************************
C
      DOUBLE PRECISION B1, B2, CM, CM2, CMM, D1, D2, DUMMY,
     .                 E, E1, E2, EMS, S, SB1, SB2, SIG, SM,
     .                 SM2, TM, TM1, TM2, TS
      B1 = B
      B2 = 1. - B1
C
C Compute slope S and second differences D1 and D2 scaled
C   by the separation between P1 and P2.
C
      S = H2 - H1
      D1 = S - HP1
      D2 = HP2 - S
C
C Test the range of SIGMA.
C
      SIG = ABS(SIGMA)
      IF (SIG .LT. 1.D-9) THEN
C
C Hermite cubic interpolation:
C
        CSHVAL = H1 + B2*(HP1 + B2*(D1 + B1*(D1 - D2)))
      ELSEIF (SIG .LE. .5) THEN
C
C 0 < SIG .LE. .5.  Use approximations designed to avoid
C   cancellation error in the hyperbolic functions.
C
        SB2 = SIG*B2
        CALL CSSNHCSH (SIG, SM,CM,CMM)
        CALL CSSNHCSH (SB2, SM2,CM2,DUMMY)
        E = SIG*SM - CMM - CMM
        CSHVAL = H1 + B2*HP1 + ((CM*SM2-SM*CM2)*(D1+D2) + SIG*
     .                     (CM*CM2-(SM+SIG)*SM2)*D1)/(SIG*E)
      ELSE
C
C SIG > .5.  Use negative exponentials in order to avoid
C   overflow.  Note that EMS = EXP(-SIG).
C
        SB1 = SIG*B1
        SB2 = SIG - SB1
        E1 = EXP(-SB1)
        E2 = EXP(-SB2)
        EMS = E1*E2
        TM = 1. - EMS
        TS = TM*TM
        TM1 = 1. - E1
        TM2 = 1. - E2
        E = TM*(SIG*(1.+EMS) - TM - TM)
        CSHVAL = H1 + B2*S + (TM*TM1*TM2*(D1+D2) + SIG*
     .                      ((E2*TM1*TM1-B1*TS)*D1 +
     .                       (E1*TM2*TM2-B2*TS)*D2))/(SIG*E)
      ENDIF
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20080612 17:17:32 nbsp Welcome to the Wiki Howdy, Kevin, and Welcome to the Wiki!   My names Evan, pleased to meet you!  Thanks for commenting on the Graduate Tutorial Services LLC entry, but it might be a better idea just to edit it and make it a more informative entry.  If you dont speak up and contribute what you know, who will?  Heck, I wrote most of the current content from Googling for information... everybody is trying to get accurate information, and youre in a position to assist, so it would be great if youd help out!  Once again, welcome to the wiki! Users/JabberWokky Evan JabberWokky Edwards


text:
program gen_be_diags_read
!------------------------------------------------------------------------
!  Purpose: Process background error statistics for display of various
!           diagnostics using WRFDA/graphics/NCL graphics utilities
!
!  Auothor: Syed RH Rizvi (MMM/NESL/NCAR)   Date: 02/01/2010
!           & Monika Krysta, (CTBTO, Vienna, Austria)
!
!  Note: Please acknowledge author/institute in work that uses this code.
!------------------------------------------------------------------------

   use da_control, only : do_normalize,filename_len,stderr,stdout,use_rf
   use da_tools_serial, only : da_get_unit
   use da_gen_be, only : da_print_be_stats_h_regional, &
      da_print_be_stats_h_global, da_print_be_stats_p, da_print_be_stats_v

   implicit none

   character*10        :: variable          ! Variable name
   character*80        :: var80c

   character*8         :: uh_method         ! Uh_method (power, scale)
   integer             :: n_smth_sl                  ! Number of smoothing for scale-length
   integer             :: cv_options                 ! WRFDA CV_OPTIONS
   character(len=filename_len)        :: filename          ! Input filename.
   integer             :: outunit           ! Output unit for diagnostics.
   integer             :: ni, nj, nk, nk_3d ! Dimensions read in.
   integer             :: bin_type          ! Type of bin to average over. !!!DALE ADD.
   integer             :: num_bins          ! Number of 3D bins.
   integer             :: num_bins2d        ! Number of 2D bins.
   integer             :: k                 ! Loop counter.
   integer             :: kdum              ! Dummy vertical index.
   integer             :: max_wavenumber    ! Smallest scale required (ni/2 - 1).
   real                :: lat_min, lat_max  ! Used if bin_type = 2 (degrees).
   real                :: binwidth_lat      ! Used if bin_type = 2 (degrees). !!!DALE ADD..
   real                :: binwidth_hgt      ! Used if bin_type = 2 (m). !!!DALE ADD..
   real                :: hgt_min, hgt_max  ! Used if bin_type = 2 (m).
   real                :: scale_length_ps   ! Scale length for scalar ps or ps_u
   logical             :: dummy

   integer, allocatable:: bin(:,:,:)        ! Bin assigned to each 3D point.
   integer, allocatable:: bin2d(:,:)        ! Bin assigned to each 2D point.
   real, allocatable   :: regcoeff_psi_chi(:)        ! psi/chi    regression cooefficient.
   real, allocatable   :: regcoeff_psi_t(:,:,:)   ! psi/t   regression cooefficient.
   real, allocatable   :: regcoeff_psi_ps(:,:)    ! psi/ps     regression cooefficient.
   real, allocatable   :: regcoeff_psi_rh(:,:,:)  ! psi/rh     regression cooefficient.
   real, allocatable   :: regcoeff_chi_u_t(:,:,:) ! chi_u/t regression coefficient
   real, allocatable   :: regcoeff_chi_u_ps(:,:)  ! chi_u/ps   regression coefficient
   real, allocatable   :: regcoeff_chi_u_rh(:,:,:)! chi_u/rh   regression coefficient
   real, allocatable   :: regcoeff_t_u_rh(:,:,:)  ! t_u/rh  regression coefficient
   real, allocatable   :: regcoeff_ps_u_rh(:,:)   ! ps_u/rh    regression coefficient

   real, allocatable   :: e_vec(:,:)        ! Domain-averaged eigenvectors.
   real, allocatable   :: e_val(:)          ! Domain-averaged eigenvalues.
   real, allocatable   :: e_vec_loc(:,:,:)  ! Latitudinally varying eigenvectors.
   real, allocatable   :: e_val_loc(:,:)    ! Latitudinally varying eigenvalues.
   real, allocatable   :: total_power(:)    ! Total Power spectrum.
   real, allocatable   :: scale_length(:)   ! Scale length for regional application.

   namelist / gen_be_diags_nl / cv_options, do_normalize, n_smth_sl, uh_method, use_rf

   integer :: iunit,namelist_unit

   ! Hardwire because of complicated incrementing of unit numbers writing a file
   ! each time
   integer, parameter :: ounit = 71

   stderr = 0
   stdout = 6

   uh_method = 'scale'
   n_smth_sl = 0
   cv_options = 5

   call da_get_unit(iunit)
   call da_get_unit(namelist_unit)

   open(unit=namelist_unit, file='gen_be_diags_nl.nl', &
        form='formatted', status='old', action='read')
   read(namelist_unit, gen_be_diags_nl)
   close(namelist_unit)

   filename = 'be.dat'
   print '("*** Unit=",i3,3X,"filename=",a40)',iunit, filename
   open (iunit, file = filename, form='unformatted')

   !----------------------------------------------------------------------------
   ! [1] Gather regression coefficients.
   !----------------------------------------------------------------------------

   ! Read the dimensions:
   read(iunit)ni, nj, nk
   nk_3d = nk

   allocate( bin(1:ni,1:nj,1:nk) )
   allocate( bin2d(1:ni,1:nj) )

   ! Read bin info:

   read(iunit)bin_type
   read(iunit)lat_min, lat_max, binwidth_lat
   read(iunit)hgt_min, hgt_max, binwidth_hgt
   read(iunit)num_bins, num_bins2d
   read(iunit)bin(1:ni,1:nj,1:nk)
   read(iunit)bin2d(1:ni,1:nj)

   if ( cv_options /= 7 ) then
      ! 1.1 Read in regression coefficients
      allocate  (regcoeff_psi_chi(1:num_bins))
      allocate  (regcoeff_psi_t(1:nk,1:nk,1:num_bins2d))
      allocate  (regcoeff_psi_ps(1:nk,1:num_bins2d))
      allocate  (regcoeff_psi_rh(1:nk,1:nk,1:num_bins2d))
      allocate  (regcoeff_chi_u_t(1:nk,1:nk,1:num_bins2d))
      allocate  (regcoeff_chi_u_ps(1:nk,1:num_bins2d))
      allocate  (regcoeff_chi_u_rh(1:nk,1:nk,1:num_bins2d))
      allocate  (regcoeff_t_u_rh(1:nk,1:nk,1:num_bins2d))
      allocate  (regcoeff_ps_u_rh(1:nk,1:num_bins2d))

      regcoeff_psi_chi = 0.
      regcoeff_psi_t   = 0.
      regcoeff_psi_ps  = 0.
      regcoeff_psi_rh  = 0.
      regcoeff_chi_u_t = 0.
      regcoeff_chi_u_ps= 0.
      regcoeff_chi_u_rh= 0.
      regcoeff_t_u_rh  = 0.
      regcoeff_ps_u_rh = 0.

      if ( cv_options == 5 ) then
         read (iunit) regcoeff_psi_chi
         read (iunit) regcoeff_psi_ps
         read (iunit) regcoeff_psi_t
      else
         do k = 1 , 9
            read (iunit) var80c
            select case( trim(adjustl(var80c)) )
            case ('regcoeff_psi_chi')
               read (iunit) regcoeff_psi_chi
            case ('regcoeff_psi_t')
               read (iunit) regcoeff_psi_t
            case ('regcoeff_psi_ps')
               read (iunit) regcoeff_psi_ps
            case ('regcoeff_psi_rh')
               read (iunit) regcoeff_psi_rh
            case ('regcoeff_chi_u_t')
               read (iunit) regcoeff_chi_u_t
            case ('regcoeff_chi_u_ps')
               read (iunit) regcoeff_chi_u_ps
            case ('regcoeff_chi_u_rh')
               read (iunit) regcoeff_chi_u_rh
            case ('regcoeff_t_u_rh')
               read (iunit) regcoeff_t_u_rh
            case ('regcoeff_ps_u_rh')
               read (iunit) regcoeff_ps_u_rh
            case default;
               write(6,fmt='(A,A)')'Error in reading regression coefficients for: ',trim(adjustl(var80c))
               stop
            end select
         end do
      end if

      outunit = ounit + 100
      call da_print_be_stats_p( outunit, ni, nj, nk, num_bins, num_bins2d, &
                    bin, bin2d, regcoeff_psi_chi,regcoeff_psi_ps,regcoeff_psi_t)
   else
      !outunit is advanced by 3 within da_print_be_stats_p; since we don't call
      !that for cv_options=7 we need to account for it
      outunit = ounit + 103
   end if

   !----------------------------------------------------------------------------
   ! [2] Gather vertical error eigenvectors, eigenvalues.
   !----------------------------------------------------------------------------

   read(iunit)variable
   read(iunit)nk, num_bins2d

   allocate( e_vec(1:nk,1:nk) )
   allocate( e_val(1:nk) )
   allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) )
   allocate( e_val_loc(1:nk,1:num_bins2d) )

   read(iunit)e_vec
   read(iunit)e_val
   read(iunit)e_vec_loc
   read(iunit)e_val_loc
   call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
                             e_vec, e_val, e_vec_loc, e_val_loc )
   read(iunit)variable
   read(iunit)nk, num_bins2d
   read(iunit)e_vec
   read(iunit)e_val
   read(iunit)e_vec_loc
   read(iunit)e_val_loc
   call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
                             e_vec, e_val, e_vec_loc, e_val_loc )

   read(iunit)variable
   read(iunit)nk, num_bins2d
   read(iunit)e_vec
   read(iunit)e_val
   read(iunit)e_vec_loc
   read(iunit)e_val_loc
   call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
                             e_vec, e_val, e_vec_loc, e_val_loc )

   read(iunit)variable
   read(iunit)nk, num_bins2d
   read(iunit)e_vec
   read(iunit)e_val
   read(iunit)e_vec_loc
   read(iunit)e_val_loc
   call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
                             e_vec, e_val, e_vec_loc, e_val_loc )

   deallocate( e_vec )
   deallocate( e_val )
   deallocate( e_vec_loc )
   deallocate( e_val_loc )

   if (uh_method /= 'power') then
      ! 2d fields: ps_u, ps:

      read(iunit)variable
      read(iunit)nk, num_bins2d

      allocate( e_vec(1:nk,1:nk) )
      allocate( e_val(1:nk) )
      allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) )
      allocate( e_val_loc(1:nk,1:num_bins2d) )

      read(iunit)e_vec
      read(iunit)e_val
      read(iunit)e_vec_loc
      read(iunit)e_val_loc
      call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
                                e_vec, e_val, e_vec_loc, e_val_loc )
      deallocate( e_vec )
      deallocate( e_val )
      deallocate( e_vec_loc )
      deallocate( e_val_loc )
   end if
   ! To assign the dimension nk for 3d fields:
   nk = nk_3d

   if (uh_method == 'power') then
       write(6,'(/a)') '[3] Gather horizontal error power spectra.'

      do k = 1, nk
         read(iunit)variable
         read(iunit)max_wavenumber, kdum
         if ( k == 1 ) allocate( total_power(0:max_wavenumber) )
         read(iunit) dummy ! to preserve namelist format
         read(iunit)total_power(:)
         call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
      end do

      do k = 1, nk
         read(iunit)variable
         read(iunit)max_wavenumber, kdum
         read(iunit) dummy ! to preserve namelist format
         read(iunit)total_power(:)
         call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
      end do

      do k = 1, nk
         read(iunit)variable
         read(iunit)max_wavenumber, kdum
         read(iunit) dummy ! to preserve namelist format
         read(iunit)total_power(:)
         call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
      end do

      do k = 1, nk
         read(iunit)variable
         read(iunit)max_wavenumber, kdum
         read(iunit) dummy ! to preserve namelist format
         read(iunit)total_power(:)
         call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
      end do

      read(iunit)variable
      read(iunit)max_wavenumber, kdum
      read(iunit) dummy ! to preserve namelist format
      read(iunit)total_power(:)
      call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )

   else if (uh_method == 'scale   ') then

      allocate (scale_length(1:nk))

      ! psi/u:
      read(iunit) variable
      read(iunit) scale_length
      call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )

      ! chi_u/v:
      read(iunit) variable
      read(iunit) scale_length
      call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )

      ! t_u/t:
      read(iunit) variable
      read(iunit) scale_length
      call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )

      ! rh:
      read(iunit) variable
      read(iunit) scale_length
      call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )

      ! ps_u/ps:
      read(iunit) variable
      read(iunit) scale_length_ps ! This is a scalar, can be ps_u or ps depending on cv_options
      write(6,'(3a,i5)')' Scale length for variable ', trim(variable), ' in unit ', outunit
      write(outunit,'(a,i4,1pe15.5)')trim(variable), 1, scale_length_ps
      outunit = outunit + 1
      write(6,*)

      deallocate (scale_length)

   else
      write(6,'(a,": uh_method=",a," not implemented.")')__FILE__,uh_method
   endif

   close(iunit)

end program gen_be_diags_read

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
function intice()
    ! Nothing because weebs are already attracted to underage anime girls
end function intice

subroutine call_the_fbi(age)
    implicit none
    integer, intent (in) :: age

    if (age < 15) then
        call execute_command_line ("call_911")
    else
       print *, "Probably just something weird happened"
    end if
end subroutine call_the_fbi

program lolicon
    implicit none
    integer :: age

    print *, "What's the girl's age?"
    read(*,*) age

    call call_the_fbi(age)
end program lolicon

max_stars_repo_name: tanafuu/lolicon

max_issues_repo_name: tanafuu/lolicon

max_forks_repo_name: tanafuu/lolicon


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                !
! THIS FILE IS THE INTERFACE OF IREGI FOR MADLOOP5 CONVENTION    !
!                                                                !
! H.-S.SHAO, 29/10/2013                                          !
!                                                                !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE INITIREGI(RECY,LOOPLIB,THR)
  USE global
  IMPLICIT NONE
  INTEGER,INTENT(IN)::LOOPLIB
  REAL(KIND(1d0)),INTENT(IN)::THR
  LOGICAL,INTENT(IN)::RECY
  !RECYCLING=.TRUE.! USE RECYCLING, THEN FOR EACH PS, CLEAR CACHING PLEASE
  NULLIFY(ibp_save)
  NULLIFY(pave_save)
  NULLIFY(shiftpaveden_save)
  NULLIFY(ibp_save2)
  NULLIFY(pave_save2)
  NULLIFY(shiftpaveden_save2)
  NULLIFY(xyzmatrices_save)
  NULLIFY(rsmatrices_save)
  RECYCLING=RECY
  CHECK=.FALSE.
  SELECT CASE(LOOPLIB)
  CASE(1)
     ! LOOPTOOLS, WE WILL USE QCDLOOP INSTEAD
     scalarlib=1
  CASE(2)
     ! AVH,i.e. ONELOOP
     scalarlib=2
  CASE(3)
     ! QCDLOOP
     scalarlib=1
  CASE DEFAULT
     ! USE QCDLOOP FOR OTHER CASES
     scalarlib=1
  END SELECT
  ZEROTHR_IREGI=THR
  ONSHELL_IREGI=.TRUE. ! always set on shell to cure PS
  ML5_CONVENTION=.TRUE. ! set ML5 CONVENTION TO MATCH THAT OF PJFRY++ OR CUTTOOLS
END SUBROUTINE INITIREGI
SUBROUTINE IMLOOP(CTMODE,IMODE,NLOOPLINE,NLOOPCOEFS,RANK,PDEN,M2L,MUR,TIRCOEFS,STABLE)
! CTMODE=1 ---> DO NOTHING
! CTMODE=2 ---> ROTATE THE LOOP PROPAGATOR DIRECTION,i.e D0...DN-1 to DN-1...D0
! CTMODE=4 ---> DO NOTHING BUT IN QUADRUPLE PRECISION (NOT DONE)
! CTMODE=5 ---> ROTATE THE LOOP PROPAGATOR DIRECTION IN QUADRUPLE PRECISION (NOT DONE)
! IMODE=0, IBP reduction
! IMODE=1, PaVe reduction
! IMODE=2, PaVe reduction with stablility improved by IBP reduction
  USE funlib
  USE ti_reduce
  USE cti_reduce
  USE global
  IMPLICIT NONE
  INTEGER,INTENT(IN)::NLOOPLINE,RANK,NLOOPCOEFS,CTMODE
  REAL(KIND(1d0)),DIMENSION(0:3,NLOOPLINE-1),INTENT(IN)::PDEN
  COMPLEX(KIND(1d0)),DIMENSION(NLOOPLINE),INTENT(IN)::M2L
  REAL(KIND(1d0)),INTENT(IN)::MUR
  COMPLEX(KIND(1d0)),DIMENSION(0:NLOOPCOEFS-1,3),INTENT(OUT)::TIRCOEFS
  LOGICAL,INTENT(OUT)::STABLE
  COMPLEX(KIND(1d0)),DIMENSION(0:NLOOPCOEFS-1,4)::TCOEFS
  !COMPLEX(KIND(1d0)),DIMENSION(0:NLOOPCOEFS-1,3)::TIRCOEFS2
  LOGICAL::TSTABLE
  REAL(KIND(1d0)),DIMENSION(NLOOPLINE)::M2LR
  INTEGER,INTENT(IN)::IMODE
  INTEGER::nlenrank,i,j
  REAL(KIND(1d0)),DIMENSION(NLOOPLINE,0:3)::PL
  !    REAL(KIND(1d0)),PARAMETER::pi=3.141592653589793d0
  REAL(KIND(1d0)),PARAMETER::eulergamma=0.5772156649015329d0
  REAL(KIND(1d0)),PARAMETER::eumlog4pi=-1.9538085820677579324d0 ! eulergamma-log(4*pi)
  REAL(KIND(1d0)),PARAMETER::euplogpi=1.7219455507509331d0 ! eulergamma+log(pi)
  ! (6 EulerGamma^2 + \[Pi]^2 - 12 EulerGamma Log[4 \[Pi]] + 6 Log[4 \[Pi]]^2)/12
  REAL(KIND(1d0)),PARAMETER::exp4eulergamma=2.7311510211049246100d0
  ! (6 EulerGamma^2 + \[Pi]^2 + 12 EulerGamma Log[\[Pi]] + 6 Log[\[Pi]]^2)/12
  REAL(KIND(1d0)),PARAMETER::exp4eulergamma2=2.3050152732995803d0
  COMPLEX(KIND(1d0)),PARAMETER::cnorm=DCMPLX(0d0,-0.1013211836423377714438794632d0) ! 1/pi**2/i
  INTEGER::nnonzeromass,ml5_sign,ipostemp
  LOGICAL::exist
  LOGICAL::comp_mass_scheme
  INTEGER::scal_save
  !    LOGICAL::INIT=.TRUE.
  !    SAVE INIT
  !    IF(INIT)THEN
  !       RECYCLING=.TRUE. ! USE RECYCLING, THEN FOR EACH PS, CLEAR CACHING PLEASE
  !       INIT=.FALSE.
  !    ENDIF
  nlenrank=0
  DO I=0,RANK
     nlenrank=nlenrank+(3+I)*(2+I)*(1+I)/6
  ENDDO
  IF(nlenrank.GT.NLOOPCOEFS)THEN
     STOP "ERROR IN IMLOOP FOR INCORRECT RANK OR NLOOPCOEFS"
  ENDIF
  nnonzeromass=0
  massref(1:MAXNLOOP_IREGI)=0d0
  comp_mass_scheme=.FALSE.
  DO i=1,NLOOPLINE
     M2LR(i)=DREAL(M2L(i))
     IF(ABS(DIMAG(M2L(i)))/MAX(ABS(M2L(i)),1D-2).GT.1D-15)THEN
        comp_mass_scheme=.TRUE.
     ENDIF
     IF(ABS(M2LR(i)).GT.ZEROTHR_IREGI)THEN
        exist=.FALSE.
        DO J=1,nnonzeromass
           IF(ABS(M2LR(I)-massref(J))/M2LR(I).LT.ZEROTHR_IREGI)THEN
              exist=.TRUE.
              EXIT
           ENDIF
        ENDDO
        IF(.NOT.exist)THEN
           nnonzeromass=nnonzeromass+1
           IF(nnonzeromass.GT.20)STOP 'ERROR:TOO MANY INDEPENDENT MASSES !!!!'
           massref(nnonzeromass)=M2LR(I)
        ENDIF
     ENDIF
     IF(i.EQ.1)THEN
        PL(i,0:3)=0d0
     ELSE
        PL(i,0:3)=PDEN(0:3,i-1)
     ENDIF
  ENDDO
  NMASS=nnonzeromass
  STABLE_IREGI=.TRUE.
  ! IMODE,NLOOPLINE,MAXRANK,NCOEFS,PDEN,M2L,MU,TICOEFS,TSTABLE
  ! IMODE=0, IBP reduction
  ! IMODE=1, PaVe reduction
  ! IMODE=2, PaVe reduction with stablility improved by IBP reduction
  ! I HARD-CODED IMODE=2 HERE
  !IMODE=2
  IF(.NOT.comp_mass_scheme)THEN
     CALL tensor_integral_reduce(IMODE,NLOOPLINE,RANK,nlenrank,&
          PL,M2LR,MUR,TCOEFS(0:nlenrank-1,1:4),TSTABLE,CTMODE)
  ELSE
     scal_save=scalarlib
     scalarlib=2
     ! complex mass scheme
     CALL comp_tensor_integral_reduce(IMODE,NLOOPLINE,RANK,nlenrank,&
          PL,M2L,MUR,TCOEFS(0:nlenrank-1,1:4),TSTABLE,CTMODE)
     scalarlib=scal_save
  ENDIF
  STABLE=TSTABLE
  ! TO MATCH THE CONVENTION OF PJFRY++
  ! MULTIPLY GAMMA(1-eps)*(Pi)^eps instead of GAMMA(1-eps)*(4*Pi)^eps
  DO i=0,nlenrank-1
     ipostemp=POS2RANK(i)
     IF(MOD(ipostemp,2).EQ.0)THEN
        ml5_sign=1
     ELSE
        ml5_sign=-1
     ENDIF
     TIRCOEFS(i,3)=TCOEFS(i,3)*cnorm*ml5_sign
     ! GAMMA(1-eps)*(4*Pi)^eps
     !TIRCOEFS2(i,2)=(TCOEFS(i,1)+TCOEFS(i,2)+TCOEFS(i,3)*eumlog4pi)*cnorm*ml5_sign
     !TIRCOEFS2(i,1)=(TCOEFS(i,4)+(TCOEFS(i,1)+TCOEFS(i,2))*eumlog4pi&
     !     +TCOEFS(i,3)*exp4eulergamma)*cnorm*ml5_sign
     ! GAMMA(1-eps)*(Pi)^eps
     TIRCOEFS(i,2)=(TCOEFS(i,1)+TCOEFS(i,2)+TCOEFS(i,3)*euplogpi)*cnorm*ml5_sign
     TIRCOEFS(i,1)=(TCOEFS(i,4)+(TCOEFS(i,1)+TCOEFS(i,2))*euplogpi&
          +TCOEFS(i,3)*exp4eulergamma2)*cnorm*ml5_sign
  ENDDO
  ! SORT THE ORDER OF COEFFICIENTS (OLD ML5 CONV) WITH THE NEW ML5 CONV
  !CALL SORT_IREGICOEFS(RANK,nlenrank,TIRCOEFS2(0:nlenrank-1,1:3),&
  !     TIRCOEFS(0:nlenrank-1,1:3))
  RETURN
END SUBROUTINE IMLOOP
SUBROUTINE IREGI_FREE_PS
  USE global
  USE binary_tree
  IMPLICIT NONE
  CALL free_ibppave_save
  CALL free_xyzmatrices_bt(xyzmatrices_save)
  CALL free_cxyzmatrices_bt(cxyzmatrices_save)
  CALL free_rsmatrices_bt(rsmatrices_save)
  CALL free_crsmatrices_bt(crsmatrices_save)
  RETURN
END SUBROUTINE IREGI_FREE_PS


text: (fn (x:Int.x=0) => x) 1;;

max_stars_repo_name: ucsc-proglang/sage

max_issues_repo_name: ucsc-proglang/sage

max_forks_repo_name: ucsc-proglang/sage


text:
!! https://github.com/flang-compiler/flang/issues/310
!!
!! POINTER ASSIGNMENT NOT DEFINING DEFERRED LENGTH PARAMETER
!!
!! In the following example, the second pointer assignment involving
!! rank remapping is not properly defining the deferred length parameter
!! of the pointer object (see F2008 5.3.14 #3).  Key to the error is
!! that the target is a deferred-length allocatable
!!
!! $ flang flang-20171122f.f90
!! $ ./a.out
!! len1=3 (expect 3)
!! len2=65538 (expect 3)

program main
  character(:), allocatable, target :: array(:)
  character(:), pointer :: ptr1(:), ptr2(:,:)
  allocate(array(1), source=['foo'])
  ptr1 => array             ! THIS WORKS
  print '(a,2(i0,a))', 'len1=', len(ptr1), ' (expect ', len(array), ')'
  ptr2(1:1,1:1) => array    ! THIS DOES NOT WORK
  print '(a,2(i0,a))', 'len2=', len(ptr2), ' (expect ', len(array), ')'
end program

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


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

subroutine set_gmpi()
    use gprec
    use omp_lib
    use gmpi
    use ghdf5_base
    use bandstru
    implicit none
#ifdef mpi_mode
    include "mpif.h"
#endif
    integer nk1,i,imap,isum
    integer local_imp(wh%num_imp)
    character*77 f_name
    logical lexist

    f_name='GOMP.h5'
    inquire(file=f_name,exist=lexist)
    if(lexist)then
        gp%ipar=1
    endif
#ifdef mpi_mode
    if(.not.lexist)then
        f_name=file_name(gp%myrank,'GMPI')
        inquire(file=f_name,exist=lexist)
        if(lexist)gp%ipar=2
    endif
#endif
    if(gp%ipar==0)then
        allocate(gp%kvec(gp%nvec,3))
        if(kpt%dim<0)then
            write(0,'(" error in set_gmpi: kpt%dim not set yet!")')
            stop
        endif
        nk1=nint(dble(kpt%dim)/gp%nprocs)
        gp%kvec(1,1)=gp%myrank
        gp%kvec(1,3)=nk1*gp%myrank
        gp%kvec(1,2)=min(nk1,kpt%dim-gp%kvec(1,3))
    else
        call gh5_open_r(f_name,f_id)
        call gh5_read(gp%nvec,'/nvec',f_id)
        allocate(gp%kvec(gp%nvec,3))
        call gh5_read(gp%kvec,gp%nvec,3,'/KVEC',f_id)
        call gh5_close(f_id)
        if(gp%ipar==1)then
#ifdef _OPENMP
            call omp_set_num_threads(gp%nprocs)
#endif
            write(gp%io,'(" omp_num_threads=",i2)')gp%nprocs
        else
#ifdef _OPENMP
            call omp_set_num_threads(1)
#endif
        endif
    endif

    call set_kpt_diml(sum(gp%kvec(:,2)))

    !< set local_imp
    local_imp=0
    isum=0
    do i=1,wh%num_imp
        imap=wh%imap_imp(i)
        if(imap==i)then
            if(gp%myrank==mod(isum,gp%nprocs).or.gp%ipar==1)then
                local_imp(i)=i
            endif
            isum=isum+1
        else
            if(local_imp(imap)>0)local_imp(i)=i
        endif
    enddo
    call wh_set_local_imp(local_imp)
    return

end subroutine set_gmpi

subroutine output_energies(io)
    use gprec
    use warehouse
    use dcstd
    use bandstru
    implicit none
    integer,intent(in)::io

    call output_energies_wh(io)
    call output_energies_dc(io)
    if(wh%ispin==2)then
        write(io,'(" spin-up band energy = ",f0.7)')bnd%eband(1)
        write(io,'(" spin-dn band energy = ",f0.7)')bnd%eband(2)
    endif
    write(io,'(" total band energy = ",f0.7)')sum(bnd%eband)
    return

end subroutine output_energies

subroutine gh5_wrt_kswt()
    use gprec
    use warehouse
    use dcstd
    use bandstru
    use ghdf5_base
    use gmpi
    implicit none

    real(8) e_gamma_dc,eband
    integer ivec,ikp,ikpl,iks,nkp,nbands,isp
    integer nemin,nemax
    complex(8),pointer :: p_kswt(:,:)

    e_gamma_dc=sum(wh%eu2)+sum(wh%et1)-wh%edcla1-sum(dc%e)

    ikpl=0
    do ivec=1,gp%nvec
        call gh5_open_w(file_name(gp%kvec(ivec,1),'KSWT'), f_id)
        if(ivec==1)then
            call gh5_write(bnd%ef,'/e_fermi',f_id)
            call gh5_write(e_gamma_dc,'/e_gamma_dc',f_id)
            eband=sum(bnd%eband)
            call gh5_write(eband,'/e_band',f_id)
            if(wh%ispin==2)then
                call gh5_write(bnd%eband(1),'/e_band_spin1',f_id)
                call gh5_write(bnd%eband(2),'/e_band_spin2',f_id)
            endif
        endif
        do iks=1,gp%kvec(ivec,2)
            ikp=gp%kvec(ivec,3)+iks
            if(gp%ipar==1)then
                ikpl=ikp
            else
                ikpl=ikpl+1
            endif
            do isp=1,bnd%nspin_in
                nbands=bnd%ne(3,ikp,isp)-bnd%ne(2,ikp,isp)+1
                p_kswt(1:nbands,1:nbands)=>bnd%hk0(1:nbands**2,sym%ie, &
                        &ikpl,isp)
                if(isp==1)then
                    call gh5_create_group('/IKP_'//trim(int_to_str(ikp)), &
                            & f_id)
                    call gh5_write(kpt%wt(ikp), &
                            &'/IKP_'//trim(int_to_str(ikp))// &
                            &"/wt",f_id)
                endif
                call gh5_write(bnd%ne(2,ikp,isp), &
                        &'/IKP_'//trim(int_to_str(ikp))// &
                        &'/nemin_spin'//trim(int_to_str(isp)),f_id)
                call gh5_write(bnd%ne(3,ikp,isp), &
                        &'/IKP_'//trim(int_to_str(ikp))// &
                        &'/nemax_spin'//trim(int_to_str(isp)),f_id)
                call gh5_write(p_kswt,nbands,nbands, &
                        &'/IKP_'//trim(int_to_str(ikp))// &
                        &'/KSWT_SPIN'//trim(int_to_str(isp)),f_id)
            enddo
        enddo
        call gh5_close(f_id)
    enddo
    nullify(p_kswt)
    return

end subroutine gh5_wrt_kswt

max_stars_repo_name: comscope/comsuite

max_issues_repo_name: comscope/Comsuite

max_forks_repo_name: comscope/Comsuite


text:
Thumper and the Booty Groovers consits of Thumper and his Booty Groovers. They were previously known as reggae ninja (hence the myspace url). Their main goal in life is to make you dance, especially if you are an attractive female. They are a bunch of guys who love Music Scene music and are out to have a good time. And they will rock your house down if only you let them.

check them out at http://www.myspace.com/reggaeninja

Sam and his sexy self


text:
module neuralnet
    use netcdf
    implicit none
    type Dense
        integer :: input_size
        integer :: output_size
        real(kind=8), allocatable :: weights(:, :)
        real(kind=8), allocatable :: bias(:)
        character(len=10) :: activation
    end type Dense

contains

    function apply_dense(layer, input)
        ! Description: Pass a set of input data through a single dense layer and nonlinear activation function
        !
        ! Args:
        ! layer (input): a single Dense object
        ! input (input): a 2D array where the rows are different examples and
        !   the columns are different model inputs
        !
        ! Returns:
        ! The output of the dense layer as a 2D array with shape (number of inputs, number of neurons)
        type(Dense), intent(in) :: layer
        real(kind=8), dimension(:, :), intent(in) :: input
        real(kind=8), dimension(size(input, 1), layer%output_size) :: apply_dense
        integer :: i, j
        apply_dense = matmul(input, layer%weights)
        do i=1, size(input, 1)
            do j=1, size(layer%bias)
                apply_dense(i, j) = apply_dense(i, j) + layer%bias(j)
            end do
        end do
        apply_dense = apply_activation(apply_dense, layer%activation)
        return
    end function apply_dense

    function apply_activation(input, activation_type)
        ! Description: Apply a nonlinear activation function to a given array of input values.
        !
        ! Args:
        ! input: A 2D array
        ! activation_type: string describing which activation is being applied. If the activation
        !       type does not match any of the available options, the linear activation is applied.
        !       Currently supported activations are:
        !           relu
        !           elu
        !           selu
        !           sigmoid
        !           tanh
        !           softmax
        !           linear
        ! Returns:
        ! activation_type: Array of the same dimensions as input with the nonlinear activation applied.
        real(kind=8), dimension(:, :), intent(in) :: input
        character(len=10), intent(in) :: activation_type
        real(kind=8), dimension(size(input, 1)) :: softmax_sum
        real(kind=8), dimension(size(input, 1), size(input, 2)) :: apply_activation
        real(kind=8), parameter :: selu_alpha = 1.6732
        real(kind=8), parameter :: selu_lambda = 1.0507
        integer :: i, j
        select case (trim(activation_type))
            case ("relu")
                where(input < 0)
                    apply_activation = 0
                elsewhere
                    apply_activation = input
                endwhere
            case ("elu")
                where(input < 0)
                    apply_activation = exp(input) - 1
                elsewhere
                    apply_activation = input
                end where
            case ("selu")
                where(input < 0)
                    apply_activation = selu_lambda * ( selu_alpha * exp(input) - selu_alpha)
                elsewhere
                    apply_activation = selu_lambda * input
                end where
            case ("sigmoid")
                apply_activation = 1.0 / (1.0 + exp(-input))
            case ("tanh")
                apply_activation = tanh(input)
            case ("softmax")
                softmax_sum = sum(exp(input), dim=2)
                do i=1, size(input, 1)
                    do j=1, size(input, 2)
                        apply_activation(i, j) = exp(input(i, j)) / softmax_sum(i)
                    end do
                end do
            case ("linear")
                apply_activation = input
            case default
                apply_activation = input
        end select
        return
    end function apply_activation

    subroutine init_neuralnet(filename, neural_net_model)
        ! init_neuralnet
        ! Description: Loads dense neural network weights from a netCDF file and builds an array of
        ! Dense types from the weights and activations.
        !
        ! Args:
        ! filename (input): Full path to the netCDF file
        ! neural_net_model (output): array of Dense layers composing a densely connected neural network
        !
        character(len=*), intent(in) :: filename
        type(Dense), allocatable, intent(out) :: neural_net_model(:)
        integer :: ncid, num_layers_id, num_layers
        integer :: layer_names_var_id, i, layer_in_dimid, layer_out_dimid
        integer :: layer_in_dim, layer_out_dim
        integer :: layer_weight_var_id
        integer :: layer_bias_var_id
        character(len=8), allocatable :: layer_names(:)
        character(len=10) :: num_layers_dim_name = "num_layers"
        character(len=11) :: layer_name_var = "layer_names"
        character (len=11) :: layer_in_dim_name
        character (len=12) :: layer_out_dim_name
        real (kind=8), allocatable :: temp_weights(:, :)
        ! Open netCDF file
        call check(nf90_open(filename, nf90_nowrite, ncid))
        ! Get the number of layers in the neural network
        call check(nf90_inq_dimid(ncid, num_layers_dim_name, num_layers_id))
        call check(nf90_inquire_dimension(ncid, num_layers_id, &
                                          num_layers_dim_name, num_layers))
        call check(nf90_inq_varid(ncid, layer_name_var, layer_names_var_id))
        allocate(layer_names(num_layers))
        call check(nf90_get_var(ncid, layer_names_var_id, layer_names))
        print *, "load neural network " // filename
        allocate(neural_net_model(num_layers))
        ! Loop through each layer and load the weights, bias term, and activation function
        do i=1, num_layers
            layer_in_dim_name = trim(layer_names(i)) // "_in"
            layer_out_dim_name = trim(layer_names(i)) // "_out"
            layer_in_dimid = -1
            ! Get layer input and output dimensions
            call check(nf90_inq_dimid(ncid, trim(layer_in_dim_name), layer_in_dimid))
            call check(nf90_inquire_dimension(ncid, layer_in_dimid, layer_in_dim_name, layer_in_dim))
            call check(nf90_inq_dimid(ncid, trim(layer_out_dim_name), layer_out_dimid))
            call check(nf90_inquire_dimension(ncid, layer_out_dimid, layer_out_dim_name, layer_out_dim))
            call check(nf90_inq_varid(ncid, trim(layer_names(i)) // "_weights", &
                                      layer_weight_var_id))
            call check(nf90_inq_varid(ncid, trim(layer_names(i)) // "_bias", &
                                      layer_bias_var_id))
            neural_net_model(i)%input_size = layer_in_dim
            neural_net_model(i)%output_size = layer_out_dim
            ! Fortran loads 2D arrays in the opposite order from Python/C, so I
            ! first load the data into a temporary array and then apply the
            ! transpose operation to copy the weights into the Dense layer
            allocate(neural_net_model(i)%weights(layer_in_dim, layer_out_dim))
            allocate(temp_weights(layer_out_dim, layer_in_dim))
            call check(nf90_get_var(ncid, layer_weight_var_id, &
                                    temp_weights))
            neural_net_model(i)%weights = transpose(temp_weights)
            deallocate(temp_weights)
            ! Load the bias weights
            allocate(neural_net_model(i)%bias(layer_out_dim))
            call check(nf90_get_var(ncid, layer_bias_var_id, &
                                    neural_net_model(i)%bias))
            ! Get the name of the activation function, which is stored as an attribute of the weights variable
            call check(nf90_get_att(ncid, layer_weight_var_id, "activation", &
                                    neural_net_model(i)%activation))
        end do
        print *, "finished loading neural network " // filename
        call check(nf90_close(ncid))

    end subroutine init_neuralnet

    subroutine neuralnet_predict(neuralnet_model, input, prediction)
        ! neuralnet_predict
        ! Description: generate prediction from neural network model for an arbitrary set of input values
        !
        ! Args:
        ! neuralnet_model (input): Array of type(Dense) objects
        ! input (input): 2D array of input values. Each row is a separate instance and each column is a model input.
        ! prediction (output): The prediction of the neural network as a 2D array of dimension (examples, outputs)

        type(Dense), intent(in) :: neuralnet_model(:)
        real(kind=8), intent(in) :: input(:, :)
        real(kind=8), allocatable, intent(out) :: prediction(:, :)
        real(kind=8), allocatable :: hidden_in(:, :)
        real(kind=8), allocatable :: hidden_out(:, :)
        integer :: i
        allocate(hidden_in(size(input, 1), size(input, 2)))
        hidden_in = input
        do i=1, size(neuralnet_model, 1)
            allocate(hidden_out(size(input, 1), &
                                neuralnet_model(i)%output_size))
            hidden_out = apply_dense(neuralnet_model(i), hidden_in)
            deallocate(hidden_in)
            allocate(hidden_in(size(hidden_out, 1), size(hidden_out, 2)))
            hidden_in = hidden_out
            deallocate(hidden_out)
        end do
        allocate(prediction(size(hidden_in, 1), size(hidden_in, 2)))
        prediction = hidden_in
        deallocate(hidden_in)
    end subroutine neuralnet_predict

    subroutine check(status)
        ! Check for netCDF errors
        integer, intent ( in) :: status
        if(status /= nf90_noerr) then
          print *, trim(nf90_strerror(status))
          stop 2
        end if
    end subroutine check

    subroutine print_2d_array(input_array)
        ! Print 2D array in pretty format
        real(kind=8), intent(in) :: input_array(:, :)
        integer :: i, j
        do i=1, size(input_array, 1)
            do j=1, size(input_array, 2)
                write(*, fmt="(1x,a,f6.3)", advance="no") "", input_array(i, j)
            end do
            write(*, *)
        end do
    end subroutine print_2d_array
end module neuralnet

max_stars_repo_name: musicpiano/mlmicrophysics

max_issues_repo_name: musicpiano/mlmicrophysics

max_forks_repo_name: musicpiano/mlmicrophysics


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
PROGRAM testLinearSolver_Multigrid
#include "UnitTest.h"
USE UnitTest
USE IntrType
USE ExceptionHandler
USE ParameterLists
USE ParallelEnv
USE VectorTypes
USE MatrixTypes
USE LinearSolverTypes
USE LinearSolverTypes_Multigrid
USE MultigridMesh

#ifdef FUTILITY_HAVE_PETSC
#include <petscversion.h>
#if (((PETSC_VERSION_MAJOR>=3) && (PETSC_VERSION_MINOR>6)) || (PETSC_VERSION_MAJOR>=4))
USE PETSCKSP
#endif
#endif

IMPLICIT NONE

TYPE(ExceptionHandlerType),TARGET :: e
TYPE(MPI_EnvType) :: mpiTestEnv
TYPE(ParamType) :: pList, optListLS, optListMat, vecPList

#ifdef FUTILITY_HAVE_PETSC
#include <petscversion.h>
#if (((PETSC_VERSION_MAJOR>=3) && (PETSC_VERSION_MINOR>=6)) || (PETSC_VERSION_MAJOR>=4))
#include <petsc/finclude/petsc.h>
#else
#include <finclude/petsc.h>
#endif
#undef IS
PetscErrorCode  :: ierr

CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)
#endif

#ifdef HAVE_MPI
CALL mpiTestEnv%init(PE_COMM_WORLD)
#endif

!> set up default parameter list
CALL optListLS%clear()
CALL optListLS%add('LinearSolverType->TPLType',LINSYS_NATIVE)
CALL optListLS%add('LinearSolverType->solverMethod',1_SNK) ! GE or BICGSTAB
CALL optListLS%add('LinearSolverType->MPI_Comm_ID',mpiTestEnv%comm)
CALL optListLS%add('LinearSolverType->numberOMP',1_SNK)
CALL optListLS%add('LinearSolverType->timerName','LinearSolver Timer')
CALL optListLS%add('LinearSolverType->matType',SPARSE)
! set parameters for matrices
CALL optListLS%add('LinearSolverType->A->MatrixType->n',-1_SNK)
CALL optListLS%add('LinearSolverType->A->MatrixType->nnz',-1_SNK)
CALL optListLS%add('LinearSolverType->A->MatrixType->isSym',.FALSE.)
! set parameters for vectors
CALL optListLS%add('LinearSolverType->x->VectorType->n',-1_SNK)
CALL optListLS%add('LinearSolverType->b->VectorType->n',-1_SNK)

! Set up vector parameter list
CALL vecPList%add('VectorType -> n',2)
CALL vecPList%add('VectorType -> MPI_Comm_ID',mpiTestEnv%comm)

!Configure exception handler for test
CALL e%setStopOnError(.FALSE.)
CALL e%setQuietMode(.TRUE.)
CALL eParams%addSurrogate(e)
CALL eLinearSolverType%addSurrogate(e)

CREATE_TEST('Test Multigrid Linear Solvers')

REGISTER_SUBTEST('testClear',testClear)
REGISTER_SUBTEST('testInit',testInit)
#ifdef FUTILITY_HAVE_PETSC
REGISTER_SUBTEST('testPreAllocPETScInterpMat',testPreAllocPETScInterpMat)
#if ((PETSC_VERSION_MAJOR<=3) && (PETSC_VERSION_MINOR<8))
REGISTER_SUBTEST('testFillInterpMats',testFillInterpMats)
REGISTER_SUBTEST('testSetupPETScMG',testSetupPETScMG)
REGISTER_SUBTEST('testIterativeSolve_Multigrid',testIterativeSolve_Multigrid)
REGISTER_SUBTEST('testSetSmoother',testSetSmoother)
REGISTER_SUBTEST('testSoftReset',testSoftReset)
REGISTER_SUBTEST('testSetInterp',testSetInterp)
#endif
#endif

FINALIZE_TEST()

CALL pList%clear()
CALL vecPList%clear()
CALL optListMat%clear()
CALL optListLS%clear()

#ifdef FUTILITY_HAVE_PETSC
CALL PetscFinalize(ierr)
#else
CALL mpiTestEnv%finalize()
#endif
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
SUBROUTINE testClear()
  TYPE(LinearSolverType_Multigrid) :: thisLS
  LOGICAL(SBK) :: bool

  !first build one by hand to test
  thisLS%isInit=.TRUE.
  thisLS%TPLType=PETSC
  thisLS%solverMethod=1
  thisLS%info=2
  thisLS%normType=2
  thisLS%maxIters=2
  thisLS%iters=2
  thisLS%relConvTol=2._SRK
  thisLS%absConvTol=2._SRK
  thisLS%residual=2._SRK
  thisLS%isDecomposed=.TRUE.
  CALL thisLS%MPIparallelEnv%init(PE_COMM_SELF)
  CALL thisLS%OMPparallelEnv%init(1)
#ifdef FUTILITY_HAVE_PETSC
  CALL KSPCreate(thisLS%MPIparallelEnv%comm,thisLS%ksp,ierr)
#endif

  ! initialize matrix A
  ALLOCATE(DenseSquareMatrixType :: thisLS%A)
  CALL pList%clear()
  CALL pList%add('MatrixType->n',2_SNK)
  CALL pList%add('MatrixType->isSym',.TRUE.)
  CALL pList%validate(pList,optListMat)
  CALL thisLS%A%init(pList) !2x2, symmetric

  ! initialize vector X
  CALL vecPList%set('VectorType -> n',2)
  ALLOCATE(RealVectorType :: thisLS%X)
  CALL thisLS%X%init(vecPList)

  ! initialize vector b
  ALLOCATE(RealVectorType :: thisLS%b)
  CALL thisLS%b%init(vecPList)

  ! initialize matrix M
  ALLOCATE(DenseSquareMatrixType :: thisLS%M)
  CALL pList%clear()
  CALL pList%add('MatrixType->n',10_SNK)
  CALL pList%add('MatrixType->isSym',.TRUE.)
  CALL thisLS%M%init(pList)

  CALL thisLS%clear()

  !check results
  bool=thisLS%isInit .OR.thisLS%solverMethod == 1                     &
      .OR. ALLOCATED(thisLS%M) .OR. ASSOCIATED(thisLS%A)              &
      .OR. ASSOCIATED(thisLS%X) .OR. thisLS%info /= 0                 &
      .OR. thisLS%normType == 2 .OR. thisLS%maxIters == 2             &
      .OR. thisLS%iters == 2 .OR. thisLS%isDecomposed                 &
      .OR. thisLS%residual == 2._SRK .OR. thisLS%relConvTol == 2._SRK &
      .OR. ALLOCATED(thisLS%PreCondType) .OR. thisLS%nLevels /= 1 &
      .OR. ALLOCATED(thisLS%level_info)                               &
      .OR. ALLOCATED(thisLS%level_info_local)                         &
#ifdef FUTILITY_HAVE_PETSC
      .OR. ASSOCIATED(thisLS%interpMats_PETSc)                        &
#endif
      .OR. thisLS%isMultigridSetup
  ASSERT(.NOT.(bool),'CALL Multigrid%clear() FAILED!')
  CALL thisLS%clear()

ENDSUBROUTINE testClear
!
!-------------------------------------------------------------------------------
SUBROUTINE testInit()
  TYPE(LinearSolverType_Multigrid) :: thisLS
  INTEGER(SIK) :: ref_level_info(2,4)
  LOGICAL(SBK) :: bool

  CALL pList%clear()
  CALL pList%add('LinearSolverType->matType',SPARSE)
  CALL pList%add('LinearSolverType->TPLType',PETSC)
  CALL pList%add('LinearSolverType->solverMethod',MULTIGRID)
  CALL pList%add('LinearSolverType->MPI_Comm_ID',PE_COMM_SELF)
  CALL pList%add('LinearSolverType->numberOMP',1_SNK)
  CALL pList%add('LinearSolverType->A->MatrixType->n',65_SNK)
  CALL pList%add('LinearSolverType->A->MatrixType->nlocal',65_SNK)
  CALL pList%add('LinearSolverType->A->MatrixType->nnz',193_SNK)
  CALL pList%add('LinearSolverType->x->VectorType->n',65_SNK)
  CALL pList%add('LinearSolverType->x->VectorType->nlocal',65_SNK)
  CALL pList%add('LinearSolverType->b->VectorType->n',65_SNK)
  CALL pList%add('LinearSolverType->b->VectorType->nlocal',65_SNK)
  CALL pList%add('LinearSolverType->PCType','PCMG')
  CALL pList%add('LinearSolverType->PCIters',-1)
  CALL pList%add('LinearSolverType->PCSetup',0)
  CALL pList%add('LinearSolverType->Multigrid->nx_local',65)
  CALL pList%add('LinearSolverType->Multigrid->ny_local',1)
  CALL pList%add('LinearSolverType->Multigrid->nz_local',1)
  CALL pList%add('LinearSolverType->Multigrid->num_eqns',1)
  CALL pList%add('LinearSolverType->Multigrid->manuallySetLevelInfo',.FALSE.)

  CALL pList%validate(pList,optListLS)
  CALL thisLS%init(pList)
  bool = (thisLS%isInit .AND. thisLS%solverMethod == MULTIGRID &
      .AND. thisLS%MPIparallelEnv%isInit() &
      .AND. thisLS%OMPparallelEnv%isInit() )
  ASSERT(bool, 'Iterative%init(...)')

  ref_level_info = RESHAPE((/1,9, &
                             1,17, &
                             1,33, &
                             1,65/),SHAPE(ref_level_info))
  ASSERT(thisLS%nLevels == 4,'Check number of multigrid levels')
  ASSERT(ALL(ref_level_info == thisLS%level_info),'Check grid sizes on each level.')

  CALL thisLS%clear()

ENDSUBROUTINE testInit
!
!-------------------------------------------------------------------------------
#ifdef FUTILITY_HAVE_PETSC
SUBROUTINE testPreAllocPETScInterpMat()
  TYPE(LinearSolverType_Multigrid) :: thisLS
  INTEGER(SIK),PARAMETER :: ref_sizes(4)=(/9,17,33,65/)
  INTEGER(SIK) :: iLevel,m,n
  PetscErrorCode :: iperr

  IF(mpiTestEnv%master) THEN
    CALL init_MultigridLS(thisLS)
    CALL preAllocInterpMatrices_1D1G(thisLS)
    DO iLevel=1,thisLS%nLevels-1
      CALL MatGetSize(thisLS%interpMats_PETSc(iLevel)%a,m,n,iperr)
      ASSERT(m == ref_sizes(iLevel+1),'Incorrect # of rows for interp. mat.')
      FINFO() 'Interpolation from',iLevel,'to',iLevel+1
      ASSERT(n == ref_sizes(iLevel),'Incorrect # of cols for interp. mat.')
      FINFO() 'Interpolation from',iLevel,'to',iLevel+1
    ENDDO
    CALL thisLS%clear()
  ELSE
    !This is needed because the wrapper assumes that if one processor has
    ! >0 asserts, then all procs have >0 asserts.  Without this, it hangs
    ASSERT(.TRUE.,'blank assert')
  ENDIF
  CALL mpiTestEnv%barrier()

ENDSUBROUTINE testPreAllocPETScInterpMat
#endif
!
!-------------------------------------------------------------------------------
#ifdef FUTILITY_HAVE_PETSC
SUBROUTINE testFillInterpMats()
  TYPE(LinearSolverType_Multigrid) :: thisLS
  TYPE(MultigridMeshStructureType),POINTER :: myMMeshes => NULL()
  TYPE(InterpWeightsStructureType),POINTER :: myWtStructure => NULL()

  INTEGER(SIK) :: iLevel,ix,ncol,nx
  INTEGER(SIK),ALLOCATABLE :: cols(:),tmpint(:)
  REAL(SRK),ALLOCATABLE :: vals(:)
  CHARACTER(LEN=2) :: tmpchar

  LOGICAL(SBK) :: boolcols,boolvals

  PetscErrorCode :: iperr

  CALL init_MultigridLS(thisLS)

  ALLOCATE(myMMeshes)
  CALL myMMeshes%init(thisLS%nLevels)
  DO iLevel=thisLS%nLevels,1,-1
    nx=thisLS%level_info(2,iLevel)
    myMMeshes%meshes(iLevel)%istt=1
    myMMeshes%meshes(iLevel)%istp=nx
    myMMeshes%meshes(iLevel)%nPointsLocal=nx
    ALLOCATE(myMMeshes%meshes(iLevel)%mmData(1:nx))
    ALLOCATE(myMMeshes%meshes(iLevel)%interpDegrees(1:nx))
    DO ix=1,nx
      IF(MOD(ix,2) == 1) THEN
        myMMeshes%meshes(iLevel)%interpDegrees(ix)=0
        ALLOCATE(myMMeshes%meshes(iLevel)%mmData(ix)%childIndices(1))
        myMMeshes%meshes(iLevel)%mmData(ix)%childIndices(1)=(ix+1)/2
      ELSE
        myMMeshes%meshes(iLevel)%interpDegrees(ix)=1
        ALLOCATE(myMMeshes%meshes(iLevel)%mmData(ix)%childIndices(2))
        myMMeshes%meshes(iLevel)%mmData(ix)%childIndices(1)=ix+1
        myMMeshes%meshes(iLevel)%mmData(ix)%childIndices(2)=ix-1
      ENDIF
    ENDDO
  ENDDO
  ALLOCATE(myWtStructure)
  CALL myWtStructure%init(myMMeshes,1)
  CALL thisLS%fillInterpMats(myMMeshes,myWtStructure)
  ALLOCATE(tmpint(thisLS%nLevels))
  tmpint(1)=GMRES
  tmpint(2:thisLS%nLevels)=SOR
  CALL pList%add('LinearSolverType->smootherMethod_list',tmpint)
  DEALLOCATE(tmpint)
  CALL thisLS%setupPETScMG(pList)

  ALLOCATE(vals(2),cols(2))
  cols=0
  DO iLevel=1,thisLS%nLevels-1
    boolcols=.TRUE.
    boolvals=.TRUE.
    nx=thisLS%level_info(2,iLevel)
    DO ix=1,nx
      vals=0.0_SRK
      CALL MatGetRow(thisLS%interpMats_PETSc(iLevel)%a,ix-1,ncol,cols,vals,iperr)
      IF(MOD(ix,2) == 1) THEN
        boolcols=(cols(1)==(ix-1)/2)
        boolvals=vals(1)==1.0_SRK
      ELSE
        boolcols=(cols(1)==ix/2-1) .AND. (cols(2)==ix/2)
        boolvals=ALL(vals(1:2)==0.5_SRK)
      ENDIF
      CALL MatRestoreRow(thisLS%interpMats_PETSc(iLevel)%a,1_SIK,ncol,cols,vals,iperr)
      IF(.NOT. boolcols .OR. .NOT. boolvals .OR. iperr /= 0 ) EXIT
    ENDDO
    WRITE(tmpchar,'(I2)') iLevel
    ASSERT(boolcols,'Check interpolation matrix col indices for grid '//tmpchar)
    ASSERT(boolvals,'Check interpolation matrix entries for grid '//tmpchar)
    ASSERT(iperr == 0,'Error obtaining matrix entries for grid '//tmpchar)

  ENDDO
  CALL thisLS%clear()
  CALL myWtStructure%clear()
  CALL myMMeshes%clear()
  DEALLOCATE(vals,cols,myWtStructure,myMMeshes)

ENDSUBROUTINE testFillInterpMats
#endif
!
!-------------------------------------------------------------------------------
#ifdef FUTILITY_HAVE_PETSC
SUBROUTINE testSetupPETScMG()
  TYPE(LinearSolverType_Multigrid) :: thisLS
  PC :: pc_temp
  KSPType :: myksptype
  PCType :: mypctype
  PetscErrorCode :: iperr

  CALL init_MultigridLS(thisLS)
  CALL preAllocInterpMatrices_1D1G(thisLS)
  CALL setupInterpMatrices_1D1G(thisLS)
  CALL thisLS%setupPETScMG(pList)
  ASSERT(thisLS%isMultigridSetup,'LS%isMultigridSetup')

  CALL KSPGetType(thisLS%ksp,myksptype,iperr)
  ASSERT(myksptype == KSPRICHARDSON,'KSP type must be richardson for MG solver in PETSc')
  CALL KSPGetPC(thisLS%ksp,pc_temp,iperr)
  CALL PCGetType(pc_temp,mypctype,iperr)
  ASSERT(mypctype == PCMG,'PC type should be Multigrid!')

  CALL thisLS%clear()
  CALL mpiTestEnv%barrier()

ENDSUBROUTINE testSetupPETScMG
#endif
!
!-------------------------------------------------------------------------------
SUBROUTINE testIterativeSolve_Multigrid()
#ifdef FUTILITY_HAVE_PETSC
  TYPE(LinearSolverType_Multigrid) :: thisLS
  REAL(SRK),ALLOCATABLE :: soln(:)
  REAL(SRK),POINTER :: x(:)
  LOGICAL(SBK) :: match
  INTEGER(SIK) :: istt,istp

  INTEGER(SIK),PARAMETER :: n=65_SNK
  INTEGER(SIK),PARAMETER :: n_2G=260_SNK

  INTEGER(SIK) :: level_info(2,4),level_info_local(2,4)

  !================ 1G Problem ============================
#ifdef HAVE_MPI
  IF(mpiTestEnv%master) THEN
#endif
    CALL init_MultigridLS(thisLS)

    ! Create solution:
    ALLOCATE(soln(n))
    soln=1.0_SRK
    soln(n/3)=2.0_SRK
    CALL setupLinearProblem_1D1G(thisLS,soln)

    CALL preAllocInterpMatrices_1D1G(thisLS)
    CALL setupInterpMatrices_1D1G(thisLS)
    CALL thisLS%setupPETScMG(pList)

    ! build x0
    ALLOCATE(x(n))
    x(1:(n-1)/2)=0.5_SRK
    x((n+1)/2:n)=1.1_SRK
    CALL thisLS%setX0(x)

    !set iterations and convergence information and build/set M
    CALL thisLS%setConv(2,1.0E-9_SRK,1.0E-9_SRK,1000,30)

    !solve it
    CALL thisLS%solve()

    SELECTTYPE(LS_x => thisLS%X); TYPE IS(PETScVectorType)
      CALL LS_x%getAll(x)
    ENDSELECT
    match=ALL(ABS(x-soln) < 1.0E-6_SRK)
    ASSERT(match, 'PETScIterative%solve() - 1G,1D Multigrid')

    DEALLOCATE(soln)
    DEALLOCATE(x)
    CALL thisLS%A%clear()
    CALL thisLS%clear()
#ifdef HAVE_MPI
  ENDIF
  CALL mpiTestEnv%barrier()
#endif
  !================ 1G Problem ============================
  !
  !
  !================ 2G,2-proc Problem ============================
#ifdef HAVE_MPI
  IF(mpiTestEnv%nproc == 2) THEN
    level_info=RESHAPE((/2,18,2,34,2,66,2,130/),(/2,4/))
    level_info_local=RESHAPE((/2,9,2,17,2,33,2,65/),(/2,4/))
    CALL init_MultigridLS(thisLS,num_eqns_in=2,nx_in=130, &
      nprocs_in=2,level_info=level_info, &
      level_info_local=level_info_local,nlevels=4)

    ! Create solution:
    ALLOCATE(soln(n_2G))
    soln=1.0_SRK
    soln(n_2G/3)=2.0_SRK
    CALL setupLinearProblem_1D2G_2proc(thisLS,soln)

    CALL preAllocInterpMatrices_1D2G_2proc(thisLS)
    CALL setupInterpMatrices_1D2G_2proc(thisLS)
    CALL thisLS%setupPETScMG(pList)

    ! build x0
    ALLOCATE(x(n_2G))
    x(1:n_2G/2)=1.1_SRK
    x(n_2G/2+1:n_2G)=0.5_SRK
    IF(mpiTestEnv%master) THEN
      istt=1
      istp=n_2G/2
    ELSE
      istt=n_2G/2+1
      istp=n_2G
    ENDIF
    SELECTTYPE(LS_x => thisLS%X); TYPE IS(PETScVectorType)
      CALL LS_x%setRange_array(istt,istp,x(istt:istp))
    ENDSELECT

    !set iterations and convergence information and build/set M
    CALL thisLS%setConv(2,1.0E-9_SRK,1.0E-9_SRK,1000,30)

    !solve it
    CALL mpiTestEnv%barrier()
    CALL thisLS%solve()

    x=0.0_SRK
    SELECTTYPE(LS_x => thisLS%X); TYPE IS(PETScVectorType)
      CALL LS_x%getRange(istt,istp,x(istt:istp))
    ENDSELECT
    match=ALL(ABS(x(istt:istp)-soln(istt:istp)) < 1.0E-6_SRK)
    ASSERT(match, 'PETScIterative%solve() - 2G,1D Multigrid')

    DEALLOCATE(soln)
    DEALLOCATE(x)
    CALL thisLS%A%clear()
    CALL thisLS%clear()
  ENDIF
#endif
  !================ 2G,2-proc Problem ============================
#endif
  !This is needed because if one or more procs have at least one assert
  ! statement, the unit test code assumes all procs have at least one
  ! assert statement.
  ASSERT(.TRUE.,'Dummy assert to make sure no procs are hanging.')

ENDSUBROUTINE testIterativeSolve_Multigrid
!
!-------------------------------------------------------------------------------
SUBROUTINE testSetSmoother
#ifdef FUTILITY_HAVE_PETSC
  TYPE(LinearSolverType_Multigrid) :: thisLS
  LOGICAL(SBK) :: tmpbool
  INTEGER(SIK) :: iLevel
  KSP :: ksp_temp
  PC :: pc_temp
  KSPType :: myksptype
  PCType :: mypctype
  PetscErrorCode :: iperr

  CALL init_MultigridLS(thisLS)
  CALL preAllocInterpMatrices_1D1G(thisLS)
  CALL setupInterpMatrices_1D1G(thisLS)
  CALL thisLS%setupPETScMG(pList)

  CALL thisLS%setSmoother(SOR,0)
  CALL PCMGGetSmoother(thisLS%pc,0,ksp_temp,iperr)
  CALL KSPGetPC(ksp_temp,pc_temp,iperr)
  CALL KSPGetType(ksp_temp,myksptype,iperr)
  CALL PCGetType(pc_temp,mypctype,iperr)
  tmpbool=(myksptype == KSPRICHARDSON) .AND. (mypctype == PCSOR)
  ASSERT(tmpbool,'Set smoother for coarsest level.')

  CALL thisLS%setSmoother(GMRES)
  tmpbool=.TRUE.
  DO iLevel=1,thisLS%nLevels-1
    CALL PCMGGetSmoother(thisLS%pc,iLevel,ksp_temp,iperr)
    CALL KSPGetPC(ksp_temp,pc_temp,iperr)
    CALL KSPGetType(ksp_temp,myksptype,iperr)
    CALL PCGetType(pc_temp,mypctype,iperr)
    tmpbool=tmpbool .AND. (myksptype == KSPGMRES)
  ENDDO
  ASSERT(tmpbool,'Set smoother for all but the coarsest level.')

  CALL thisLS%setSmoother(SOR,1)
  CALL PCMGGetSmoother(thisLS%pc,1,ksp_temp,iperr)
  CALL KSPGetPC(ksp_temp,pc_temp,iperr)
  CALL KSPGetType(ksp_temp,myksptype,iperr)
  CALL PCGetType(pc_temp,mypctype,iperr)
  tmpbool=(myksptype == KSPRICHARDSON) .AND. (mypctype == PCSOR)
  DO iLevel=2,thisLS%nLevels-1
    CALL PCMGGetSmoother(thisLS%pc,iLevel,ksp_temp,iperr)
    CALL KSPGetPC(ksp_temp,pc_temp,iperr)
    CALL KSPGetType(ksp_temp,myksptype,iperr)
    CALL PCGetType(pc_temp,mypctype,iperr)
    tmpbool=tmpbool .AND. (myksptype == KSPGMRES)
  ENDDO
  ASSERT(tmpbool,'Set smoother level 1 only.')

  CALL thisLS%clear()
#endif

ENDSUBROUTINE testSetSmoother
!
!-------------------------------------------------------------------------------
SUBROUTINE testSoftReset
#ifdef FUTILITY_HAVE_PETSC
  TYPE(LinearSolverType_Multigrid) :: thisLS
  REAL(SRK),ALLOCATABLE :: soln(:)
  REAL(SRK),POINTER :: x(:)
  TYPE(ParamType) :: matPList
  INTEGER(SIK),ALLOCATABLE :: dnnz(:),onnz(:)
  INTEGER(SIK) :: i
  LOGICAL(SBK) :: match
  INTEGER(SIK),PARAMETER :: n=65_SNK
  !================ 1G Problem ============================
#ifdef HAVE_MPI
  IF(mpiTestEnv%master) THEN
#endif
    CALL init_MultigridLS(thisLS)

    ! Create solution:
    ALLOCATE(soln(n))
    soln=1.0_SRK
    soln(n/3)=2.0_SRK
    CALL setupLinearProblem_1D1G(thisLS,soln)

    CALL preAllocInterpMatrices_1D1G(thisLS)
    CALL setupInterpMatrices_1D1G(thisLS)
    CALL thisLS%setupPETScMG(pList)

    ! build x0
    ALLOCATE(x(n))
    x(1:(n-1)/2)=0.5_SRK
    x((n+1)/2:n)=1.1_SRK
    CALL thisLS%setX0(x)

    !set iterations and convergence information and build/set M
    CALL thisLS%setConv(2,1.0E-9_SRK,1.0E-9_SRK,1000,30)

    !solve it
    CALL thisLS%solve()

    !prepare matrix A for reinit
    CALL thisLS%A%clear()
    DEALLOCATE(thisLS%A)

    !set up parameter list
    CALL matPList%clear()
    CALL matPList%add('MatrixType->engine',VM_PETSC)
    CALL matPList%add('MatrixType->matType',SPARSE)
    CALL matPList%add('MatrixType->MPI_Comm_ID',PE_COMM_SELF)
    CALL matPList%add('MatrixType->n',n)
    CALL matPList%add('MatrixType->nlocal',n)
    ALLOCATE(dnnz(n),onnz(n))
    onnz=0
    dnnz=1
    CALL matPList%add('MatrixType->dnnz',dnnz)
    CALL matPList%add('MatrixType->onnz',onnz)
    DEALLOCATE(dnnz,onnz)

    !allocate and initialize matrix (A)
    CALL MatrixFactory(thisLS%A,matPList)
    !set A to the identity matrix:
    DO i=1,n
      CALL thisLS%A%set(i,i,1.0_SRK)
      CALL thisLS%b%set(i,1.0_SRK)
    ENDDO
    !associate A with the KSP
    CALL thisLS%updatedA()
    !clean up
    CALL matPList%clear()

    ! build x0
    DEALLOCATE(x)
    ALLOCATE(x(n))
    x(1:(n-1)/2)=0.5_SRK
    x((n+1)/2:n)=1.1_SRK
    CALL thisLS%setX0(x)

    !solve new system
    CALL thisLS%solve()
    SELECTTYPE(LS_x => thisLS%X); TYPE IS(PETScVectorType)
      CALL LS_x%getAll(x)
    ENDSELECT
    match=ALL(ABS(x-1.0_SRK) < 1.0E-6_SRK)
    ASSERT(match, 'testSoftReset()')

    DEALLOCATE(soln)
    DEALLOCATE(x)
    CALL thisLS%clear()

#ifdef HAVE_MPI
  ELSE
    !This is needed or it will hang:
    ASSERT(.TRUE., 'dummy assert for proc 2')
  ENDIF
  CALL mpiTestEnv%barrier()
#endif
#endif

ENDSUBROUTINE testSoftReset
!
!-------------------------------------------------------------------------------
! This subroutine tests solving two linear systems with the same interp. matrices
SUBROUTINE testSetInterp
#ifdef FUTILITY_HAVE_PETSC
  TYPE(LinearSolverType_Multigrid) :: thisLS,thisLS2
  REAL(SRK),ALLOCATABLE :: soln(:),soln2(:)
  REAL(SRK),POINTER :: x(:)
  LOGICAL(SBK) :: match

  INTEGER(SIK),PARAMETER :: n=65_SNK

#ifdef HAVE_MPI
  IF(mpiTestEnv%master) THEN
#endif
    CALL init_MultigridLS(thisLS)
    CALL init_MultigridLS(thisLS2)

    ! Create solution:
    ALLOCATE(soln(n))
    soln=1.0_SRK
    soln(n/3)=2.0_SRK
    CALL setupLinearProblem_1D1G(thisLS,soln)
    ! Create solution:
    ALLOCATE(soln2(n))
    soln2=1.0_SRK
    CALL setupLinearProblem_1D1G(thisLS2,soln2)

    CALL preAllocInterpMatrices_1D1G(thisLS)
    CALL setupInterpMatrices_1D1G(thisLS)
    ! Reuse interp matrices:
    CALL thisLS2%setInterp(thisLS%interpMats_PETSc)

    CALL thisLS%setupPETScMG(pList)
    CALL thisLS2%setupPETScMG(pList)

    !set iterations and convergence information and build/set M
    CALL thisLS%setConv(2,1.0E-9_SRK,1.0E-9_SRK,1000,30)
    CALL thisLS2%setConv(2,1.0E-9_SRK,1.0E-9_SRK,1000,30)

    ! build x0
    ALLOCATE(x(n))
    x(1:(n-1)/2)=0.5_SRK
    x((n+1)/2:n)=1.1_SRK
    CALL thisLS%setX0(x)
    CALL thisLS2%setX0(x)
    !solve it
    CALL thisLS%solve()
    CALL thisLS2%solve()
    SELECTTYPE(LS_x => thisLS%X); TYPE IS(PETScVectorType)
      CALL LS_x%getAll(x)
    ENDSELECT
    match=ALL(ABS(x-soln) < 1.0E-6_SRK)
    SELECTTYPE(LS_x => thisLS2%X); TYPE IS(PETScVectorType)
      CALL LS_x%getAll(x)
    ENDSELECT
    match=match .AND. ALL(ABS(x-soln2) < 1.0E-6_SRK)
    ASSERT(match, 'testSetInterp()')

    DEALLOCATE(soln2)
    DEALLOCATE(soln)
    DEALLOCATE(x)
    CALL thisLS%A%clear()
    CALL thisLS%clear()
    CALL thisLS2%A%clear()
    CALL thisLS2%clear()
#ifdef HAVE_MPI
  ELSE
    !This is needed or it will hang:
    ASSERT(.TRUE., 'dummy assert for proc 2')
  ENDIF
  CALL mpiTestEnv%barrier()
#endif
#endif
ENDSUBROUTINE testSetInterp
!
!-------------------------------------------------------------------------------
SUBROUTINE init_MultigridLS(thisLS,num_eqns_in,nx_in,ny_in,nz_in,nprocs_in, &
      level_info,level_info_local,nlevels)
  TYPE(LinearSolverType_Multigrid),INTENT(INOUT) :: thisLS
  INTEGER(SIK),INTENT(IN),OPTIONAL :: nx_in,ny_in,nz_in,num_eqns_in,nprocs_in
  INTEGER(SIK),INTENT(IN),OPTIONAL :: level_info(:,:),level_info_local(:,:)
  INTEGER(SIK),INTENT(IN),OPTIONAL :: nlevels
  INTEGER(SIK) :: num_eqns,nx,ny,nz,nprocs,n,nlocal
  INTEGER(SIK),ALLOCATABLE :: dnnz(:),onnz(:)

  nprocs=1
  num_eqns=1
  nx=65
  ny=1
  nz=1
  IF(PRESENT(num_eqns_in)) num_eqns=num_eqns_in
  IF(PRESENT(nx_in)) nx=nx_in
  IF(PRESENT(ny_in)) ny=ny_in
  IF(PRESENT(nz_in)) nz=nz_in
  IF(PRESENT(nprocs_in)) nprocs=nprocs_in

  n=num_eqns*nx*ny*nz
  nlocal=n/nprocs

  !The PETSC sparse matrix type
  ! initialize linear system
  CALL pList%clear()
  CALL pList%add('LinearSolverType->TPLType',PETSC)
  CALL pList%add('LinearSolverType->solverMethod',MULTIGRID)
  IF(nprocs > 1) THEN
    CALL pList%add('LinearSolverType->MPI_Comm_ID',mpiTestEnv%comm)
  ELSE
    CALL pList%add('LinearSolverType->MPI_Comm_ID',PE_COMM_SELF)
  ENDIF
  CALL pList%add('LinearSolverType->numberOMP',1_SNK)
  CALL pList%add('LinearSolverType->timerName','testTimer')
  CALL pList%add('LinearSolverType->matType',SPARSE)
  CALL pList%add('LinearSolverType->A->MatrixType->n',n)
  CALL pList%add('LinearSolverType->A->MatrixType->nlocal',nlocal)
  ALLOCATE(dnnz(nlocal),onnz(nlocal))
  onnz=0
  dnnz(1:num_eqns)=1+num_eqns
  dnnz(num_eqns+1:nlocal-num_eqns)=2+num_eqns
  dnnz(nlocal-num_eqns+1:nlocal)=1+num_eqns
  IF(nprocs > 1) THEN
    IF(mpiTestEnv%master) THEN
      onnz(nlocal-num_eqns+1:nlocal)=1
    ELSE
      onnz(1:num_eqns)=1
    ENDIF
  ENDIF
  CALL pList%add('LinearSolverType->A->MatrixType->dnnz',dnnz)
  CALL pList%add('LinearSolverType->A->MatrixType->onnz',onnz)
  CALL pList%add('LinearSolverType->A->MatrixType->nnz',SUM(dnnz)+SUM(onnz))
  DEALLOCATE(dnnz,onnz)
  CALL pList%add('LinearSolverType->x->VectorType->n',n)
  CALL pList%add('LinearSolverType->x->VectorType->nlocal',nlocal)
  CALL pList%add('LinearSolverType->b->VectorType->n',n)
  CALL pList%add('LinearSolverType->b->VectorType->nlocal',nlocal)
  CALL pList%validate(pList,optListLS)

  !Geometry dimensions for multigrid:
  CALL pList%add('LinearSolverType->Multigrid->nx_local',nx/nprocs)
  CALL pList%add('LinearSolverType->Multigrid->ny_local',ny)
  CALL pList%add('LinearSolverType->Multigrid->nz_local',nz)
  CALL pList%add('LinearSolverType->Multigrid->num_eqns',num_eqns)
  IF(PRESENT(nlevels)) &
      CALL pList%add('LinearSolverType->Multigrid->nlevels',nlevels)
  IF(PRESENT(level_info)) &
      CALL pList%add('LinearSolverType->Multigrid->level_info',level_info)
  IF(PRESENT(level_info_local)) &
      CALL pList%add('LinearSolverType->Multigrid->level_info_local', &
      level_info_local)
  IF(.NOT.(PRESENT(nlevels) .AND. PRESENT(level_info) .AND. &
      PRESENT(level_info_local))) &
      CALL pList%add('LinearSolverType->Multigrid->manuallySetLevelInfo',.FALSE.)

  CALL thisLS%init(pList)

ENDSUBROUTINE init_MultigridLS
!
!-------------------------------------------------------------------------------
SUBROUTINE preAllocInterpMatrices_1D1G(thisLS)
  TYPE(LinearSolverType_Multigrid),INTENT(INOUT) :: thisLS
#ifdef FUTILITY_HAVE_PETSC
  INTEGER(SIK) :: iLevel, inx, nx, nx_old
  INTEGER(SIK),ALLOCATABLE :: dnnz(:),onnz(:)

  IF(.NOT.ASSOCIATED(thisLS%interpMats_PETSc)) THEN
    ALLOCATE(thisLS%interpMats_PETSc(thisLS%nLevels-1))
    thisLS%isOwnerOfInterpMats=.TRUE.
  ENDIF
  nx=thisLS%level_info(2,thisLS%nLevels)
  DO iLevel=thisLS%nLevels-1,1,-1
    !Create the interpolation operator:
    nx_old=nx
    nx=thisLS%level_info(2,iLevel)

    ALLOCATE(dnnz(nx_old),onnz(nx_old))
    onnz=0
    DO inx=1,nx_old
      IF(XOR(MOD(inx,2)==1, (inx>nx/2 .AND. MOD(nx,2) == 0))) THEN
        dnnz(inx)=1
      ELSE
        dnnz(inx)=2
      ENDIF
    ENDDO

    CALL thisLS%preAllocPETScInterpMat(iLevel,dnnz,onnz)
    DEALLOCATE(dnnz,onnz)

  ENDDO
#endif
ENDSUBROUTINE preAllocInterpMatrices_1D1G
!
!-------------------------------------------------------------------------------
SUBROUTINE setupInterpMatrices_1D1G(thisLS)
  TYPE(LinearSolverType_Multigrid),INTENT(INOUT) :: thisLS
#ifdef FUTILITY_HAVE_PETSC
  INTEGER(SIK) :: iLevel, inx, nx, nx_old

  nx=thisLS%level_info(2,thisLS%nLevels)
  DO iLevel=thisLS%nLevels-1,1,-1
    !Create the interpolation operator:
    nx_old=nx
    nx=thisLS%level_info(2,iLevel)

    DO inx=1,nx_old
      IF(XOR(MOD(inx,2)==1, (inx>nx/2 .AND. MOD(nx,2) == 0))) THEN
        CALL thisLS%interpMats_PETSc(iLevel)%set(inx,(inx+1)/2,1.0_SRK)
      ELSE
        CALL thisLS%interpMats_PETSc(iLevel)%set(inx,inx/2,0.5_SRK)
        CALL thisLS%interpMats_PETSc(iLevel)%set(inx,inx/2+1,0.5_SRK)
      ENDIF
    ENDDO

  ENDDO
#endif
ENDSUBROUTINE setupInterpMatrices_1D1G
!
!-------------------------------------------------------------------------------
SUBROUTINE preAllocInterpMatrices_1D2G_2proc(thisLS)
  TYPE(LinearSolverType_Multigrid),INTENT(INOUT) :: thisLS
#ifdef FUTILITY_HAVE_PETSC
  INTEGER(SIK) :: iLevel,inx,nx,nx_old,ieqn,row
  INTEGER(SIK),ALLOCATABLE :: dnnz(:),onnz(:)

  IF(.NOT.ASSOCIATED(thisLS%interpMats_PETSc)) THEN
    ALLOCATE(thisLS%interpMats_PETSc(thisLS%nLevels-1))
    thisLS%isOwnerOfInterpMats=.TRUE.
  ENDIF
  nx=thisLS%level_info_local(2,thisLS%nLevels)
  DO iLevel=thisLS%nLevels-1,1,-1
    !Create the interpolation operator:
    nx_old=nx
    nx=thisLS%level_info_local(2,iLevel)

    ALLOCATE(dnnz(nx_old*2),onnz(nx_old*2))
    onnz=0
    DO inx=1,nx_old
      DO ieqn=1,2
        row=(inx-1)*2+ieqn
        IF(XOR(MOD(inx,2)==1, (inx>nx/2 .AND. MOD(nx,2) == 0))) THEN
          dnnz(row)=1
        ELSE
          dnnz(row)=2
        ENDIF
      ENDDO
    ENDDO

    CALL thisLS%preAllocPETScInterpMat(iLevel,dnnz,onnz)
    DEALLOCATE(dnnz,onnz)

  ENDDO
#endif
ENDSUBROUTINE preAllocInterpMatrices_1D2G_2proc
!
!-------------------------------------------------------------------------------
SUBROUTINE setupInterpMatrices_1D2G_2proc(thisLS)
  TYPE(LinearSolverType_Multigrid),INTENT(INOUT) :: thisLS
#ifdef FUTILITY_HAVE_PETSC
  INTEGER(SIK) :: iLevel,inx,nx,nx_old,ieqn,row
  INTEGER(SIK) :: inx_c,col
  INTEGER(SIK) :: offset(2)

  nx=thisLS%level_info_local(2,thisLS%nLevels)
  DO iLevel=thisLS%nLevels-1,1,-1
    !Create the interpolation operator:
    nx_old=nx
    nx=thisLS%level_info_local(2,iLevel)
    IF(mpiTestEnv%master) THEN
      offset=0
    ELSE
      offset(1)=nx_old*2
      offset(2)=nx*2
    ENDIF

    DO inx=1,nx_old
      DO ieqn=1,2
        row=(inx-1)*2+ieqn
        IF(XOR(MOD(inx,2)==1, (inx>nx/2 .AND. MOD(nx,2) == 0))) THEN
          col=(inx-1)+ieqn
          CALL thisLS%interpMats_PETSc(iLevel)%set(row+offset(1),col+offset(2),1.0_SRK)
        ELSE
          inx_c=inx/2
          col=(inx_c-1)*2+ieqn
          CALL thisLS%interpMats_PETSc(iLevel)% &
              set(row+offset(1),col+offset(2),0.5_SRK)
          col=inx_c*2+ieqn
          CALL thisLS%interpMats_PETSc(iLevel)% &
              set(row+offset(1),col+offset(2),0.5_SRK)
        ENDIF
      ENDDO
    ENDDO

  ENDDO
#endif
ENDSUBROUTINE setupInterpMatrices_1D2G_2proc
!
!-------------------------------------------------------------------------------
SUBROUTINE setupLinearProblem_1D1G(thisLS,soln)
  TYPE(LinearSolverType_Multigrid),INTENT(INOUT) :: thisLS
  REAL(SRK),INTENT(IN) :: soln(:)
#ifdef FUTILITY_HAVE_PETSC
  REAL(SRK),ALLOCATABLE :: b(:)
  REAL(SRK),ALLOCATABLE :: A_temp(:,:)
  INTEGER(SIK) :: n
  INTEGER(SIK) :: i

  n=SIZE(soln)

  !A is a tridiagonal system with -1 on the offdiagonals, and
  !  2.5 on the diagonals.
  ALLOCATE(A_temp(n,n))
  A_temp=0.0_SRK
  SELECTTYPE(A => thisLS%A); TYPE IS(PETScMatrixType)
    DO i=1,n
      IF(i > 1) THEN
        CALL A%set(i,i-1,-1.0_SRK)
        A_temp(i,i-1)=-1.0_SRK
      ENDIF
      CALL A%set(i,i,2.5_SRK)
      A_temp(i,i)=2.5_SRK
      IF(i < n) THEN
        CALL A%set(i,i+1,-1.0_SRK)
        A_temp(i,i+1)=-1.0_SRK
      ENDIF
    ENDDO
  ENDSELECT

  ALLOCATE(b(n))
  b=MATMUL(A_temp,soln)
  DEALLOCATE(A_temp)

  SELECTTYPE(LS_b => thisLS%b); TYPE IS(PETScVectorType)
    CALL LS_b%setAll_array(b)
  ENDSELECT
  DEALLOCATE(b)
#endif

ENDSUBROUTINE setupLinearProblem_1D1G
!
!-------------------------------------------------------------------------------
SUBROUTINE setupLinearProblem_1D2G_2proc(thisLS,soln)
  TYPE(LinearSolverType_Multigrid),INTENT(INOUT) :: thisLS
  REAL(SRK),INTENT(IN) :: soln(:)
#ifdef FUTILITY_HAVE_PETSC
#ifdef HAVE_MPI
  REAL(SRK),ALLOCATABLE :: b(:)
  REAL(SRK),ALLOCATABLE :: A_temp(:,:)
  INTEGER(SIK) :: n,nstart,nend
  INTEGER(SIK) :: i,grp
  INTEGER(SIK) :: row,col

  n=SIZE(soln)/2

  !2G homogeneous diffusion problem:
  ALLOCATE(A_temp(2*n,2*n))
  A_temp=0.0_SRK
  SELECTTYPE(A => thisLS%A); TYPE IS(PETScMatrixType)
    IF(mpiTestEnv%master) THEN
      nstart=1
      nend=n/2
    ELSE
      nstart=n/2+1
      nend=n
    ENDIF
    DO i=nstart,nend
      DO grp=1,2
        row=(i-1)*2+grp
        IF(i > nstart) THEN
          col=(i-2)*2+grp
          CALL A%set(row,col,-1.0_SRK)
          A_temp(row,col)=-1.0_SRK
        ENDIF
        IF(i < nend) THEN
          col=i*2+grp
          CALL A%set(row,col,-1.0_SRK)
          A_temp(row,col)=-1.0_SRK
        ENDIF
      ENDDO
      row=(i-1)*2+1
      CALL A%set(row,row,2.1_SRK)
      A_temp(row,row)=2.1_SRK
      CALL A%set(row+1,row+1,2.2_SRK)
      A_temp(row+1,row+1)=2.2_SRK
      !Upscatter:
      CALL A%set(row,row+1,-0.01_SRK)
      A_temp(row,row+1)=-0.01_SRK
      !Downscatter:
      CALL A%set(row+1,row,-0.5_SRK)
      A_temp(row+1,row)=-0.5_SRK
    ENDDO
  ENDSELECT

  CALL mpiTestEnv%allReduce(4*n*n,A_temp)

  ALLOCATE(b(2*n))
  b=MATMUL(A_temp,soln)
  DEALLOCATE(A_temp)

  SELECTTYPE(LS_b => thisLS%b); TYPE IS(PETScVectorType)
    CALL LS_b%setRange_array((nstart-1)*2+1,(nend-1)*2+2, &
        b((nstart-1)*2+1:(nend-1)*2+2))
  ENDSELECT
  DEALLOCATE(b)
#endif
#endif

ENDSUBROUTINE setupLinearProblem_1D2G_2proc

ENDPROGRAM testLinearSolver_Multigrid

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: wgurecky/Futility


SUBROUTINE UpdateGhostLayerX(var,S,Nx,Ny,DiffStencils,alpha,GhostGridX,GhostGridY)
! By Allan P. Engsig-Karup.
USE Precision
USE DataTypes
IMPLICIT NONE
TYPE (Diff_def) :: DiffStencils
INTEGER :: Nx, Ny, rank, alpha, i,j, GhostGridX, GhostGridY
REAL(KIND=long), DIMENSION(Nx,Ny) :: var, S
REAL(KIND=long), DIMENSION(2*alpha+1) :: Stencil
rank = 2*alpha+1
! two boundaries
i=1
DO j = 1+GhostGridY, Ny-GhostGridY
    Stencil = DiffStencils%StencilX(i+GhostGridX,1:rank,1)
    var(i,j) = (S(i+GhostGridX,j)-DOT_PRODUCT(Stencil(2:rank),var(2:rank,j) )) / Stencil(1)
END DO
i = Nx
DO j = 1+GhostGridY, Ny-GhostGridY
    Stencil = DiffStencils%StencilX(i-GhostGridX,1:rank,1)
    var(i,j) = (S(i-GhostGridX,j)-DOT_PRODUCT(Stencil(1:rank-1),var(i-rank+1:i-1,j) )) / Stencil(rank)
END DO
END SUBROUTINE UpdateGhostLayerX


text:
../tb/Driver.sv
../tb/Interface.sv
../tb/Monitor.sv
../tb/Generator.sv
../tb/Scoreboard.sv
../tb/Environment.sv

../tb/Testbench.sv
../tb/clock_reset.sv

../tc/iic_1/testcase.sv

../../../digital/rtl/chip_top.v
../../../digital/rtl/chip.v
../../../digital/rtl/clk_gen.v

../../../digital/rtl/bus/bus.v
../../../digital/rtl/bus/bus_addr_dec.v
../../../digital/rtl/bus/bus_arbiter.v
../../../digital/rtl/bus/bus_master_mux.v
../../../digital/rtl/bus/bus_slave_mux.v

../../../digital/rtl/cpu/alu.v
../../../digital/rtl/cpu/bus_if.v
../../../digital/rtl/cpu/cpu.v
../../../digital/rtl/cpu/ctrl.v
../../../digital/rtl/cpu/decoder.v
../../../digital/rtl/cpu/ex_reg.v
../../../digital/rtl/cpu/ex_stage.v
../../../digital/rtl/cpu/gpr.v
../../../digital/rtl/cpu/id_reg.v
../../../digital/rtl/cpu/id_stage.v
../../../digital/rtl/cpu/if_reg.v
../../../digital/rtl/cpu/if_stage.v
../../../digital/rtl/cpu/mem_ctrl.v
../../../digital/rtl/cpu/mem_reg.v
../../../digital/rtl/cpu/mem_stage.v
../../../digital/rtl/cpu/spm.v
../../../digital/rtl/cpu/symmetric_mem_core.v

../../../digital/rtl/gpio/gpio.v
../../../digital/rtl/mem/rom.v
../../../digital/rtl/timer/timer.v

../../../digital/rtl/uart/uart.v
../../../digital/rtl/uart/uart_ctrl.v
../../../digital/rtl/uart/uart_tx.v
../../../digital/rtl/uart/uart_rx.v


text:
!***********************************************************************
      subroutine sys_mkdir (dir, lendir, ierr)
!
!   This routine makes a sub-dir under the current working directory.
!   lendir is the length of character string dir;
!   ierr will be zero if successful, otherwise non-zero;
!   machine is an optional parameter specifying the name of the system
!
!   Xinghong He  98-08-21
!
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:58   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!
!   function 'system' and execute_command_line commented out
!cjb                                             Jacek Bieron 2018 June 18
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
!     use mpi_C
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      implicit none
      character(len=*), intent(in):: dir
      integer, intent(in):: lendir
      integer, intent(out):: ierr
      character(len=255) mkdir_error_message
!      character(len=*), parameter::iam = 'SUN'  ! soly for machine
!      character(len=len_trim (iam)), optional, intent(out):: machine

!     integer system

      mkdir_error_message = " no message passed from the shell "
      ierr = -999

!cjb
!     print*, dir(1:lendir)
!     print*, dir(1:lendir)
!
      ierr = system ('mkdir -p -m 775 ' // dir(1:lendir))
!cjb  use EXECUTE_COMMAND_LINE if function 'system' is not supported
!     call execute_command_line ('mkdir -p -m 775 ' // dir(1:lendir), &
!                         exitstat=ierr, cmdmsg=mkdir_error_message)
      if (ierr.ne.0) then
        print*, ' sys_mkdir error message = ', trim(mkdir_error_message)
!               ' myid = ', myid
      endif
      return
      end subroutine sys_mkdir

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
! ------------------------------------------------------------
!     This program illustrates the following points:
!     (1)  The exponential trap:
!          That is, A**B**C is equal to A**(B**C) rather
!          than (A**B)**C.
!     (2)  The integer division trap:
!          That is, 4/6 is ZERO in Fortran rather than
!          a real number 0.666666
!          Function REAL() is used to illustrate the
!          differences.
!     (3)  The string truncation trap:
!          What if the length assigned to a CHARACTER
!          is shorter than the length of the string you
!          expect the identifier to have?  The third part
!          shows you the effect.
! ------------------------------------------------------------

! gfortran -c fortran_traps.f90
! gfortran fortran_traps.o -o fortran_traps

PROGRAM fortran_traps
   IMPLICIT NONE

   INTEGER, PARAMETER          :: A = 2, B = 2, H = 3
   INTEGER, PARAMETER          :: O = 4, P = 6
   CHARACTER(LEN=5), PARAMETER :: M = 'Smith', N = 'TEXAS'
   CHARACTER(LEN=4), PARAMETER :: X = 'Smith'
   CHARACTER(LEN=6), PARAMETER :: Y = 'TEXAS'

!  The exponential trap

   WRITE(*,*)   "First, the exponential trap:"
   WRITE(*,*)   A, ' ** ', B, ' ** ', H, ' = ', A**B**H
   WRITE(*,*)   '( ', A, ' ** ', B, ' ) **', H, ' = ', (A**B)**H
   WRITE(*,*)   A, ' ** ( ', B, ' ** ', H, ' ) = ', A**(B**H)
   WRITE(*,*)

!  The integer division trap.  Intrinsic function REAL() converts
!  an integer to a real number

   WRITE(*,*)   "Second, the integer division trap:"
   WRITE(*,*)
   WRITE(*,*)   O, ' / ', P, ' = ', O/P
   WRITE(*,*)   'REAL( ', O, ' ) / ', P, ' = ', REAL(O)/P
   WRITE(*,*)   O, ' / REAL( ', P, ' ) = ', O/REAL(P)
   WRITE(*,*)

!  The string truncation trap

   WRITE(*,*)   "Third, the string truncation trap:"
   WRITE(*,*)   'IS ', M, ' STILL IN ', N, '?'
   WRITE(*,*)   'IS ', X, ' STILL IN ', Y, '?'

END PROGRAM fortran_traps

max_stars_repo_name: NickPepper/Fortran

max_issues_repo_name: NickPepper/Fortran

max_forks_repo_name: NickPepper/Fortran