text:
	SUBROUTINE GTXSY ( sys, itype, isym, ijust, ixoff, iyoff,
     +			   rotn, x, y, cchar, iret )
C************************************************************************
C* GTXSY								*
C*									*
C* This subroutine plots a special text and symbol to any coordinate    *
C* system.  The special text is centered on the given reference point   *
C* (X,Y).  The text is drawn using attributes defined by GSTEXT, and    *
C* the surrounding box is drawn with attributes defined by GSLINE.      *
C* Depending upon special text type, the box may be filled.             *
C*                                                                      *
C* GTXSY ( SYS, ITYPE, ISYM, IJUST, IXOFF, IYOFF, ROTN, X, Y,           *
C*	   CCHAR, IRET )						*
C*									*
C* Input parameters:							*
C*	SYS		CHAR*		Coordinate system		*
C*                                        'S' = screen coordinates      *
C*					  'D' = device coordinates	*
C*					  'N' = normalized coordinates	*
C*					  'V' = view coordinates	*
C*					  'P' = plot coordinates	*
C*					  'M' = map coordinates		*
C*					  'G' = grid coordinates	*
C*	ITYPE		INTEGER		Special Text Type		*
C*					    1 = right justified		*
C*                                          2 = high pressure box       *
C*                                          3 = pressure box            *
C*                                                  bounded, unfilled   *
C*                                          4 = pressure box            *
C*                                                  bounded, filled     *
C*                                          5 = pressure box,           *
C*                                                  unbounded, filled   *
C*                                          6 = freezing level symbol   *
C*					    7 = turbulence symbol	*
C*					    8 = cloud level		*
C*					    9 = high level turbulence	*
C*					   10 = underline		*
C*					   11 = underline, fill box	*
C*					   12 = midlevel icing		*
C*					   13 = overline		*
C*					   14 = overline, fill box	*
C*					   15 = "Big Box" for mid-level	*
C*	ISYM		INTEGER		Symbol type if itype of 7	*
C*					    one or two digits in range  *
C*					     0-8. 			*
C*	IJUST		INTEGER		Justification (-1, 0, 1) 	*
C*	IXOFF		INTEGER		X Offset			*
C*	IYOFF		INTEGER		Y Offset			*
C*	ROTN		REAL		Rotation			*
C*	X		REAL		X coordinate / latitude		*
C*	Y		REAL		Y coordinate / longitude	*
C*	CCHAR		CHAR*		Text string to plot		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* E. Safford/GSC	 4/97	Initial coding -- copied from GTEXTC    *
C* E. Safford/GSC        5/97   Changed sym to isym                     *
C* C. Lin/EAI	 	 6/97	Modified for adding 'S' to sysup    	*
C* S. Schotz/NCEP        7/97   Update documentation for S coordinate   *
C* E. Safford/GSC        7/97   Returned isym as parameter              *
C* E. Safford/GSC        7/97   Added ijust, irotn, ioff		*
C* S. Jacobs/NCEP	 7/98	Added check for N-relative rotation	*
C* S. Jacobs/NCEP	 7/98	Redefined N-relative angles		*
C* S. Jacobs/NCEP	 7/98	Added Grid coord to N-rel check		*
C* S. Jacobs/NCEP	 9/98	Fixed rotation variable name typo	*
C* I. Durham/GSC	10/98	Changed GAZDRM to GP_AZDR		*
C* M. Li/SAIC		11/01	Added type 12 for midlevel icing	*
C* T. Lee/SAIC		 8/02	Added text to prolog			*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'DEVSET.CMN'
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'DEVACT.CMN'
C*
	CHARACTER*(*)	cchar, sys
C
	CHARACTER	sysuc*1, sysout*1
C
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = NORMAL
C
C*	Check that device has been set.
C
	IF  ( ddev .eq. ' ' ) THEN
	    iret = NDVICE
	    RETURN
	END IF
C
C*	Find coordinate system for clipping.
C
	CALL ST_LCUC ( sys, sysuc, ier )
	isys = INDEX ( sysup, sysuc )
	IF  ( isys .eq. 0 ) THEN
	    iret = NOCORD
	    RETURN
	  ELSE
	    iwndw = isys - 2
	    IF ( iwndw .gt. 3 ) iwndw = 3
	    IF ( iwndw .le. 0 ) iwndw = 1
	END IF
C
C*	Translate to device coordinates.
C
	IF  ( ddev .eq. 'VG' )  THEN
	    sysout = 'M'
	  ELSE
	    sysout = 'D'
	END IF
C
C*	Check for North-relative rotation.
C
	IF  ( ( mrrotn .eq. 2 )  .and.
     +	      ( ( sysuc .eq. 'M' ) .or. ( sysuc .eq. 'G' ) ) .and.
     +	      ( ddev .ne. 'VG' ) )  THEN
C
C*	    First convert from Grid to Map coordinates, if needed.
C
	    IF  ( sysuc .eq. 'G' )  THEN
		CALL GTRANS ( sysuc, 'M', 1, x, y, xout, yout, ier )
		sysuc = 'M'
	      ELSE
		xout = x
		yout = y
	    END IF
C
C*	    Apply the North relative axis rotation.
C
	    CALL GP_AZDR ( rotn, xout, yout, rangl, ier )
C
C*	    Compute the angle as if it were for Screen relative text.
C
	    rotat = rotn + ( rotn - rangl )
	  ELSE
	    rotat = rotn
	    xout  = x
	    yout  = y
	END IF
C
C*	Transform to the output coordinate system.
C
	CALL GTRANS ( sysuc, sysout, 1, xout, yout, xp, yp, ier )
C
C*	Check for missing location.
C
	IF ( ( ERMISS ( xp ) ) .or. ( ERMISS ( yp ) ) ) RETURN
C
C*	Reset color if necessary.
C
	IF ( lcolr .ne. mcolr ) CALL DSCOLR ( lcolr, jcm, ier )
C
C*	Send the text string to the device driver.
C
	CALL ST_LSTR ( cchar, lenc, ier )
	CALL DTXSY  ( iwndw, itype, isym, ijust, ixoff, iyoff, rotat,
     +		      xp, yp, cchar, lenc, iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! PR fortran/26246
! { dg-options "-fdump-tree-original" }
! { dg-do compile }

module pr26246_1
  implicit none
  contains
    function foo(string)
    character(*), intent(in) :: string
    character(len=len(string)+2) :: foo
    if (index(trim(string), '"').ne.0) then
      foo = "'" // trim(string) // "'"
    end if
  end function foo
end module pr26246_1

! { dg-final { scan-tree-dump-times "static int" 0 "original" } }
! { dg-final { cleanup-tree-dump "original" } }
! { dg-final { cleanup-modules "pr26246_1" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      subroutine dumpdatex(ifile)
C
C #####################################################################
C
C     PURPOSE -
C
C        None
C
C     INPUT ARGUMENTS -
C
C        None
C
C     OUTPUT ARGUMENTS -
C
C        None
C
C     CHANGE HISTORY -
C
C        $Log: dumpdatex.f,v $
C        Revision 2.00  2007/11/05 19:45:53  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.0   08 Feb 2000 09:42:40   dcg
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.12   Tue Mar 31 13:08:58 1998   dcg
CPVCS    fix format statements for newer version of datex standard
CPVCS
CPVCS       Rev 1.11   Wed Jul 16 10:23:44 1997   dcg
CPVCS    declare logmess
CPVCS
CPVCS       Rev 1.10   Mon Apr 14 16:44:02 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.9   Tue Mar 05 12:49:14 1996   dcg
CPVCS    remove icn1, int1
CPVCS
CPVCS       Rev 1.8   Fri Feb 23 16:36:02 1996   dcg
CPVCS    remove explicit references to uic,vic,wic,pic,ric,eic
CPVCS
CPVCS       Rev 1.7   11/07/95 17:16:18   dcg
CPVCS    change flag to 2 in mmgetblk calls
CPVCS
CPVCS       Rev 1.6   06/20/95 15:43:10   dcg
CPVCS    remove character literals from arguments list to hgetprt
CPVCS
CPVCS       Rev 1.5   05/26/95 13:13:24   het
CPVCS    Replace subroutine parameter list with subroutine calles.
CPVCS
CPVCS       Rev 1.4   05/01/95 08:35:54   het
CPVCS    Modifiy all the storage block calles for long names
CPVCS
CPVCS       Rev 1.3   03/03/95 12:09:20   dcg
CPVCS     Remove hard-wired 2dmesh mesh object in decimate, dump routines
CPVCS
CPVCS       Rev 1.2   01/04/95 22:02:06   llt
CPVCS    unicos changes (made by het)
CPVCS
CPVCS       Rev 1.1   12/19/94 08:20:44   het
CPVCS    Add the "comdict.h" include file.
CPVCS
CPVCS
CPVCS       Rev 1.0   11/13/94 11:44:18   pvcs
CPVCS    Orginal Version
C
C ######################################################################
C
      implicit none
      include 'geom_lg.h'
c
C ######################################################################
C
      character ifile*(*)
      pointer (ipimt1, imt1)
      pointer (ipitp1, itp1)
      pointer (ipicr1, icr1)
      pointer (ipisn1, isn1)
      pointer (ipicn1, icn1)
      pointer (ipxic, xic)
      pointer (ipyic, yic)
      pointer (ipzic, zic)
      pointer (ipvels,vels)
      pointer (ippres, pres)
      pointer (ipdens, dens)
      pointer (ipener, ener)
      pointer (ipitetclr, itetclr)
      pointer (ipitet, itet)
      pointer (ipjtet, jtet)
      pointer (ipout,out)
      integer nsd, nen, nef
      integer imt1(*), itp1(*),
     *        icr1(*), isn1(*), icn1(*)
      real*8 xic(*), yic(*), zic(*)
      real*8 vels(3,*)
      real*8 dens(*), pres(*)
      real*8 ener(*),rout,out(*)
      integer itetclr(*)
      integer itet(*), jtet(*)
C
      character*32 isubname, iprtname, jfile, cmo, cvelnm,
     *        cpresnm, cdensnm, cenernm,geom_name

      character*132 logmess
      logical iscontact
      integer icscode,ir,iclr,i1,icontact,index,ishape,i,it,
     *  nregions,iere,ierd,ierp,ier,itype,lout,iout,ierror,
     *  iunit,nnodes,length,icmotype,nelements,mbndry,
     *  nsdtopo,nsdgeom,icharlnf,lenfile,ilen,ityp,ierr
C
      pointer (iplocation, location)
      integer location(10000000)
C
      integer iconn(100)
C
      integer jshape(6)
      data jshape / 3, 4, 4, 5, 6, 6 /
      integer x3d_to_simul_hex(8), simul_to_x3d_hex(8)
C                             1  2  3  4  5  6  7  8
      data x3d_to_simul_hex / 1, 2, 4, 5, 3, 6, 8, 7 /
      data simul_to_x3d_hex / 1, 2, 5, 3, 4, 6, 8, 7 /

C
      lenfile=icharlnf(ifile)
      jfile=ifile(1:lenfile) // '.geo'
      iunit=-1
      call hassign(iunit,jfile,ierror)
      if (iunit.lt.0 .or. ierror.lt.0) then
        call x3d_error(isubname,'hassign bad file unit')
        goto 9999
      endif

C
C get current mesh info
C
      isubname='dumpdatex'
C
      call cmo_get_name(cmo,ierror)
C
      call cmo_get_info('nnodes',cmo,
     *                  nnodes,length,icmotype,ierror)
      call cmo_get_info('nelements',cmo,
     *                  nelements,length,icmotype,ierror)
      call cmo_get_info('mbndry',cmo,
     *                  mbndry,length,icmotype,ierror)
      call cmo_get_info('ndimensions_topo',cmo,
     *                  nsdtopo,length,icmotype,ierror)
      call cmo_get_info('ndimensions_geom',cmo,
     *                  nsdgeom,length,icmotype,ierror)
      call cmo_get_info('nodes_per_element',cmo,
     *                  nen,length,icmotype,ierror)
      call cmo_get_info('faces_per_element',cmo,
     *                  nef,length,icmotype,ierror)
      call cmo_get_info('imt1',cmo,ipimt1,ilen,ityp,ierr)
      call cmo_get_info('itp1',cmo,ipitp1,ilen,ityp,ierr)
      call cmo_get_info('icr1',cmo,ipicr1,ilen,ityp,ierr)
      call cmo_get_info('isn1',cmo,ipisn1,ilen,ityp,ierr)
      call cmo_get_info('icn1',cmo,ipicn1,ilen,ityp,ierr)
      if( ierr.ne.0) then
         write(logmess,'("No contact attribute in dumpdatex")')
         call writloga('default',0,logmess,0,ierr)
         iscontact=.false.
      else
         iscontact=.true.
      endif
      call cmo_get_info('xic',cmo,ipxic,ilen,ityp,ierr)
      call cmo_get_info('yic',cmo,ipyic,ilen,ityp,ierr)
      call cmo_get_info('zic',cmo,ipzic,ilen,ityp,ierr)
      call cmo_get_attinfo('presname',cmo,iout,rout,cpresnm,
     *   ipout,lout,itype,ier)
      call cmo_get_info(cpresnm,cmo,ippres,ilen,ityp,ierp)
      call cmo_get_attinfo('densname',cmo,iout,rout,cdensnm,
     *   ipout,lout,itype,ier)
      call cmo_get_info(cdensnm,cmo,ipdens,ilen,ityp,ierd)
      call cmo_get_attinfo('enername',cmo,iout,rout,cenernm,
     *   ipout,lout,itype,ier)
      call cmo_get_info(cenernm,cmo,ipener,ilen,ityp,iere)
      call cmo_get_info('itetclr',cmo,
     *                        ipitetclr,ilen,ityp,ierr)
      call cmo_get_info('itet',cmo,ipitet,ilen,ityp,ierr)
      call cmo_get_info('jtet',cmo,ipjtet,ilen,ityp,ierr)
C
      nsd=nsdtopo
C
      nregions=0
      do it=1,nelements
         nregions=max(nregions,itetclr(it))
      enddo
      write(iunit,"('DATEX2.1')")
      write(iunit,"('Geometry')")
      write(iunit,"(a9)") '"LaGriT1.0"'
      write(iunit,"(a28)") '"SIMUL file written by LaGriT"'
      write(iunit,*) nsd
      write(iunit,*) nnodes
      write(iunit,*) nelements
      write(iunit,*) nregions
      write(iunit,"('Node')")
      write(iunit,*)
      if(nsd.eq.3) then
         do i=1,nnodes
            write(iunit,9000) xic(i),yic(i),zic(i)
9000        format(3(1pe15.7))
         enddo
      elseif(nsd.eq.2) then
         do i=1,nnodes
            write(iunit,9000) xic(i),yic(i)
         enddo
      endif
      if(nsd.eq.2.and.nen.eq.3.and.nef.eq.3) ishape=1
      if(nsd.eq.3.and.nen.eq.3.and.nef.eq.3) ishape=1
      if(nsd.eq.2.and.nen.eq.4.and.nef.eq.4) ishape=2
      if(nsd.eq.3.and.nen.eq.4.and.nef.eq.4) ishape=3
      if(nsd.eq.3.and.nen.eq.5.and.nef.eq.5) ishape=4
      if(nsd.eq.3.and.nen.eq.6.and.nef.eq.5) ishape=5
      if(nsd.eq.3.and.nen.eq.8.and.nef.eq.6) ishape=6
      write(iunit,"('Element')")
      write(iunit,*)
      do it=1,nelements
         index=nen*(it-1)
         if(ishape.eq.6) then
            do i=1,nen
               iconn(i)=itet(index+x3d_to_simul_hex(i))-1
            enddo
         else
            do i=1,nen
               iconn(i)=itet(index+i)-1
            enddo
         endif
         icontact=0
         do i=1,nen
            i1=itet(index+i)
            if(iscontact)icontact=max(icontact,icn1(i1))
         enddo
         icontact=icontact
         iclr=itetclr(it)-1
         write(iunit,*) ishape,icontact,iclr
         write(iunit,*) (iconn(i),i=1,nen)
      enddo
      write(iunit,"('Region')")
      write(iunit,*)
      call cmo_get_attinfo('geom_name',cmo,iout,rout,geom_name,
     *   ipout,lout,itype,ierr)
      call mmfindbk('cmregs',geom_name,ipcmregs,lout,ierr)
      do ir=1,nregions
         iprtname=cmregs(ir)
      enddo
      length=nnodes
      call mmgetblk("location",isubname,iplocation,length,2,icscode)
      write(iunit,"('Location')")
      do i=1,nnodes
         if (iscontact) then
            location(i)=-2
            if(itp1(i).eq.0) then
               location(i)=-2
            elseif(itp1(i).eq.2.or.itp1(i).eq.12) then
               location(i)=-1
            elseif(itp1(i).ge.10.and.itp1(i).le.19) then
               location(i)=icn1(i)
            endif

            if(icn1(i).eq.0) then
               location(i)=-2
            else
               location(i)=icn1(i)
            endif
         else
            location(i)=0
         endif
      enddo
      write(iunit,"(20i5)") (location(i),i=1,nnodes)
      close(iunit)
      iunit=-1
      lenfile=icharlnf(ifile)
      jfile=ifile(1:lenfile) // '.dat'

      call hassign(iunit,jfile,ierror)
      if (iunit.lt.0 .or. ierror.lt.0) then
        call x3d_error(isubname,'hassign bad file unit')
        goto 9999
      endif

      write(iunit,"('DATEX2.1')")
      write(iunit,"('Simulation')")
      write(iunit,"(a8)") '"X3D1.0"'
      write(iunit,"(a27)") '"SIMUL file written by X3D"'
      write(iunit,*) 4
      write(iunit,"('DopingConcentration scalar node')")
      write(iunit,"('BoronConcentration scalar node')")
      write(iunit,"('TotalConcentration scalar node')")
      write(iunit,"('Velocity vector node')")
      if(ierd.eq.0) then
      write(iunit,"('Dataset')")
      do i=1,nnodes
         write(iunit,9000) dens(i)
      enddo
      write(iunit,"('0')")
      endif
      if(ierp.eq.0) then
      write(iunit,"('Dataset')")
      do i=1,nnodes
         write(iunit,9000) pres(i)
      enddo
      write(iunit,"('0')")
      endif
      if(iere.eq.0) then
      write(iunit,"('Dataset')")
      do i=1,nnodes
         write(iunit,9000) ener(i)
      enddo
      write(iunit,"('0')")
      endif
      call cmo_get_attinfo('velname',cmo,iout,rout,cvelnm,
     *   ipout,lout,itype,ier)
      call cmo_get_info(cvelnm,cmo,ipvels,length,itype,ier)
      if(ier.eq.0) then
      write(iunit,"('Dataset')")
      do i=1,nnodes
         write(iunit,9000) vels(1,i),vels(2,i),vels(3,i)
      enddo
      write(iunit,"('0')")
      endif
      close(iunit)
      goto 9999
 9999 continue
      call mmrelprt(isubname,icscode)
      return
      end

ext: f

max_forks_repo_name: daniellivingston/LaGriT


text:
cc ------------ dpmjet3.4 - authors: S.Roesler, R.Engel, J.Ranft -------
cc -------- phojet1.12-40 - authors: S.Roesler, R.Engel, J.Ranft -------
cc                                                      - oct'13 -------
cc ----------- pythia-6.4 - authors: Torbjorn Sjostrand, Lund'10 -------
cc ---------------------------------------------------------------------
cc                                  converted for use with FLUKA -------
cc                                                      - oct'13 -------

C...PYFACT
C...Multiplies histogram contents by factor.

      SUBROUTINE PYFACT(ID,F)

C...Double precision declaration.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
C...Commonblock.
      include 'inc/pybins'

C...Find initial address in memory. Multiply all contents bins.
      IF(ID.LE.0.OR.ID.GT.IHIST(1)) CALL PYERRM(28,
     &'(PYFACT:) not allowed histogram number')
      IS=INDX(ID)
      IF(IS.EQ.0) CALL PYERRM(28,
     &'(PYFACT:) scaling unbooked histogram')
      DO 100 IX=IS+6,IS+8+NINT(BIN(IS+1))
        BIN(IX)=F*BIN(IX)
  100 CONTINUE

      RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
      SUBROUTINE BLOCARF( ZS,N,LAG,NS0,KMAX,ZMEAN,SUM,AIC,C,B,A,SD,NP,
     *                    NE,SXX )
C

C
cc      PROGRAM  BLOCAR
C.......................................................................
C.....PLANNED BY H.AKAIKE...............................................
C.....DESIGNED BY H.AKAIKE AND G.KITAGAWA...............................
C.....PROGRAMMED BY G.KITAGAWA AND F.TADA...............................
C.....ADDRESS: THE INSTITUTE OF STATISTICAL MATHEMATICS, 4-6-7 MINAMI-AZ
C..............MINATO-KU, TOKYO 106, JAPAN..............................
C.....DATE OF THE LATEST REVISION:  MAR. 6,1979.........................
C.......................................................................
C.....THIS PROGRAM WAS ORIGINALLY PUBLISHED IN "TIMSAC-78", BY H.AKAIKE,
C.....G.KITAGAWA, E.ARAHATA AND F.TADA, COMPUTER SCIENCE MONOGRAPHS, NO.
C.....THE INSTITUTE OF STATISTICAL MATHEMATICS, TOKYO, 1979.............
C.......................................................................
C     TIMSAC 78.3.2.
C     _                  ___                __
C     BAYESIAN METHOD OF LOCALLY STATIONARY AR MODEL FITTING; SCALAR CAS
C
C     THIS PROGRAM LOCALLY FITS AUTOREGRESSIVE MODELS TO NON-STATIONARY
C     SERIES BY A BAYESIAN PROCEDURE.  POWER SPECTRA FOR STATIONARY SPAN
C     ARE GRAPHICALLY PRINTED OUT.  (THIS PROGRAM IS TENTATIVE.)
C
C     INPUTS REQUIRED:
C             MT:       INPUT DEVICE FOR ORIGINAL DATA (MT=5 : CARD READ
C             LAG:      UPPER LIMIT OF THE ORDER OF AR MODEL, MUST BE LE
C                       OR EQUAL TO 50.
C             NS:       LENGTH OF BASIC LOCAL SPAN
C             KSW:      =0  CONSTANT VECTOR IS NOT INCLUDED AS A REGRESS
C                       =1  CONSTANT VECTOR IS INCLUDED AS THE FIRST REG
C
C               -- THE FOLLOWING INPUTS ARE REQUESTED BY SUBROUTINE REDA
C             TITLE:    SPECIFICATION OF DATA
C             N:        DATA LENGTH, MUST BE LESS THAN OR EQUAL TO 10000
C             DFORM:    INPUT DATA SPECIFICATION STATEMENT.
C                       -- EXAMPLE  --     (8F10.5)
C             (Z(I),I=1,N):  ORIGINAL DATA
C               --------------------------------------------------------
C
C
C     AT EACH STAGE OF MODELLING OF LOCAL AR MODEL, A TWO-STEP BAYESIAN
C     PROCEDURE IS APPLIED
C        1.   AVERAGING OF THE MODELS WITH DIFFERENT ORDERS FITTED TO TH
C             NEWLY OBTAINED DATA (BY SUBROUTINE ARBAYS).
C        2.   AVERAGING OF THE MODELS FITTED TO THE PRESENT AND PRECEDIN
C             (AS FAR AS KMAX(<21)) SPANS.
C     AS TO THE AVERAGING 1, SEE THE COMMENTS OF PROGRAM UNIBAR.
C
C     AIC OF THE MODEL FITTED TO THE NEW SPAN IS DEFINED BY
C
C          AIC = NS * LOG( SD ) + 2 * EK,
C     WHERE
C          NS:     LENGTH OF NEW DATA
C          SD:     INNOVATION VARIANCE
C          EK:     EQUIVALENT NUMBER OF PARAMETERS, DEFINED AS THE SUM O
C                  SQUARES OF THE BAYESIAN WEIGHTS
C
C     AIC OF THE MODELS FITTED TO THE PRECEDING SPANS ARE DEFINED BY
C
C          AIC(J+1) = NS * LOG( SD(J) ) + 2     (J=1,KC)
C     WHERE
C          SD(J):  PREDICTION ERROR VARIANCE BY THE MODEL FITTED TO J
C                  PERIODS FORMER SPAN.
C       ---------------------------------------------------------------
C       REFERENCES:
C          G.KITAGAWA AND H.AKAIKE(1978), "A PROCEDURE FOR THE MODELING
C          OF NON-STATIONARY TIME SERIES.",  ANN. INST. STATIST. MATH.,
C          30,B,351-363.
C
C          H.AKAIKE(1978), "A BAYESIAN EXTENSION OF THE MINIMUM AIC
C          PROCEDURE OF AUTOREGRESSIVE MODEL FITTING.",  RESEARCH MEMO.
C          NO. 126, THE INSTITUTE OF STATISTICAL MATHEMATICS; TOKYO.
C       ---------------------------------------------------------------
C       THE FOLLOWING SUBROUTINES ARE DIRECTLY CALLED BY THIS PROGRAM:
C             REDATA
C             NONSTB
C             PRINTA
C             NRASPE
C       ---------------------------------------------------------------
C
cc      !DEC$ ATTRIBUTES DLLEXPORT :: BLOCARF
C
CC      IMPLICIT REAL * 8 ( A-H , O-Y )
cxx      IMPLICIT REAL * 8 ( A-H , O-Z )
cc      REAL * 4  TITLE(20) , TTL(13)
cc      DIMENSION Z(10000)
cc      DIMENSION  X(200,51) , D(200) , A(50)
cxx      DIMENSION  Z(N), ZS(N)
cxx      DIMENSION  X(NS0,LAG+1), A(LAG,KMAX)
cxx      DIMENSION  AIC(KMAX,KMAX), C(KMAX,KMAX)
cxx      DIMENSION  B(LAG,KMAX), SD(KMAX)
cxx      DIMENSION  NP(KMAX), NE(KMAX)
      INTEGER :: N, LAG, NS0, KMAX, NP(KMAX), NE(KMAX)
      REAL(8) :: ZS(N), ZMEAN, SUM, AIC(KMAX,KMAX), C(KMAX,KMAX),
     1           B(LAG,KMAX), A(LAG,KMAX), SD(KMAX), SXX(121,KMAX)
cc      DATA  TTL / 4H  CU,4HRREN,4HT MO,4HDEL ,4H(AVE,4HRAGE,4H BY ,4HTHE
cc     1 ,4HBAYE,4HSIAN,4H WEI,4HGHTS,4H)    /
C
cxx      DIMENSION    SXX(121,KMAX)
cxx      DIMENSION    F(LAG,KMAX)
      REAL(8) :: Z(N), X(NS0,LAG+1), F(LAG,KMAX), BB
C
      EXTERNAL  SETX1
C
cc      CHARACTER(100) IFLNAM,OFLNAM
cc      CALL FLNAM2( IFLNAM,OFLNAM,NFL )
cc      IF (NFL.EQ.0) GO TO 999
cc      IF (NFL.EQ.2) THEN
cc         OPEN( 6,FILE=OFLNAM,ERR=900,IOSTAT=IVAR )
cc      ELSE
cc         CALL SETWND
cc      END IF
C
C
      KSW = 0
cc      KMAX = 20
cc      MJ1 = 200
      MJ1 = NS0
      ISW = 0
C
CC      READ( 5,1 )     MT
cc      MT = 5
cc      OPEN( MT,FILE=IFLNAM,ERR=910,IOSTAT=IVAR,STATUS='OLD' )
cc      READ( 5,1 )     LAG , NS
      NS = NS0
C
cc      WRITE( 6,3 )
cc      IF( KSW .EQ. 1 )     WRITE( 6,5 )
cc      IF( KSW .NE. 1 )     WRITE( 6,4 )
cc      WRITE( 6,6 )
cc      WRITE( 6,2 )   LAG , NS , MT
C
C          ---------------------
C          ORIGINAL DATA LOADING
C          ---------------------
C
cc      CALL  REDATA( Z,N,MT,TITLE )
cc      SUBROUTINE  REDATA( X,N,MT,TITLE )
      CALL  REDATA( ZS,Z,N,ZMEAN,SUM )
cc      CLOSE( MT )
C
C     ---  INITIAL SETTING  ---
      LAG1 = LAG+1
      N0 = 0
      K = LAG + KSW
      K3 = K*3
C
      NR = 0
      KC = 0
   10 CONTINUE
      NR = NR+1
C
C          -----------------------------------
C          LOCALLY STATIONARY AR-MODEL FITTING
C          -----------------------------------
C
cc      CALL  NONSTB( SETX1,Z,X,D,LAG,N0,NS,KMAX,KSW,ISW,TITLE,MJ1,A,SD )
cx      CALL  NONSTB( SETX1,Z,X,LAG,N0,NS,KMAX,KSW,ISW,MJ1,KC,F,
      CALL  NONSTB( SETX1,Z,N,X,LAG,N0,NS,KMAX,KSW,ISW,MJ1,KC,F,
     *              AIC(1,NR),C(1,NR),B(1,NR),A(1,NR),SD(NR) )
C
cc      NP = N0 + LAG + 1
cc      NE = NP+NS-1
cc      CALL  PRINTA( A,SD,K,TTL,13,TITLE,NP,NE )
      NP(NR) = N0 + LAG + 1
      NE(NR) = NP(NR)+NS-1
C
C     ---  SPECTRUM DISPLAY  ---
cc      CALL  NRASPE( SD,A,B,K,0,120,TITLE )
      CALL  NRASPE( SD(NR),A(1,NR),BB,K,0,120,SXX(1,NR) )
C
C
C     ---  PREPARATION FOR THE NEXT STEP  ---
      N0 = N0 + NS
      IF(N0+NS+LAG1.GT.N)      NS = N-N0-LAG1
      IF(N-N0-NS-LAG1.LT.K3)   NS = N-N0-LAG1
      IF( N0+LAG1 .LT. N )     GO TO 10
C
C
cc      GO TO 999
C
cc  900 CONTINUE
cc      WRITE(6,600) IVAR,OFLNAM
cc      GO TO 999
C
cc  910 CONTINUE
cc      IF (NFL.EQ.2) CLOSE( 6 )
cc#ifdef __linux__
ccC	reopen #6 as stdout
cc      IF (NFL.EQ.2) OPEN(6, FILE='/dev/fd/1')
cc#endif
ccC /* __linux__ */
cc      WRITE(6,610) IVAR,IFLNAM
cc      GO TO 999
C
cxx  600 FORMAT(/,' !!! Output_Data_File OPEN ERROR ',I8,//5X,100A)
cxx  610 FORMAT(/,' !!! Input_Data_File OPEN ERROR ',I8,//5X,100A)
C
      RETURN
cxx    1 FORMAT( 16I5 )
cxx    2 FORMAT( ///1H ,'  FITTING UP TO THE ORDER  K =',I3,'  IS TRIED',/,
cxx     1'   BASIC LOCAL SPAN  NS =',I4,/,'   ORIGINAL DATA INPUT DEVICE  M
cxx     2T =',I3 )
cxx    3 FORMAT( //' PROGRAM TIMSAC 78.3.2',/,'   BAYESIAN METHOD OF LOCALL
cxx     *Y STATIONARY AR MODEL FITTING;   SCALAR CASE',//,'   < BASIC AUTOR
cxx     *EGRESSIVE MODEL >' )
cxx    4 FORMAT( 1H ,10X,'Z(I) = A(1)*Z(I-1) + A(2)*Z(I-2) + ... + A(M)*Z(I
cxx     1-M) + E(I)' )
cxx    5 FORMAT( 1H ,10X,'Z(I) = A(1) + A(2)*Z(I-1) + ... + A(M+1)*Z(I-M) +
cxx     1 E(I)' )
cxx    6 FORMAT( 1H ,2X,'WHERE',/,11X,'M:     ORDER OF THE MODEL',/,11X,'E(
cxx     1I):  GAUSSIAN WHITE NOISE WITH MEAN 0  AND  VARIANCE SD(M).' )
      END
cc      SUBROUTINE  NONSTB( SETX,Z,X,D,LAG,N0,NS,KMAX,KSW,ISW,TITLE,MJ1,A,
cc     1SD )
cx      SUBROUTINE  NONSTB( SETX,Z,X,LAG,N0,NS,KMAX1,KSW,ISW,MJ1,KC,F,
      SUBROUTINE  NONSTB( SETX,Z,N,X,LAG,N0,NS,KMAX1,KSW,ISW,MJ1,KC,F,
     1AIC,C,B,A,SD )
C
C     BAYESIAN TYPE NON-STATIONARY AUTOREGRESSIVE MODEL FITTING PROCEDUR
C
C     THIS SUBROUTINE FIRST FITS AN AUTOREGRESSIVE MODEL TO THE NEWLY OB
C     DATA SPAN AND THEN PRODUCES A BAYESIAN MODEL BY AVERAGING THE MODE
C     FITTED TO THE PRESENT AND PRECEDING SPANS.
C
C     BAYESIAN WEIGHT OF THE MODEL FITTED TO I PERIODS FORMER SPAN IS
C     DEFINED BY
C          C(I)  =  CONST * P(I) / (I+1)
C     WHERE
C          CONST  =  NORMALIZING CONSTANT
C          P(I)  =  EXP( -0.5*AIC(I) )
C          AIC(I)  =  AIC WITH RESPECT TO THE PRESENT DATA.
C       ----------------------------------------------------------------
C       THE FOLLOWING SUBROUTINES ARE DIRECTLY CALLED BY THIS SUBROUTINE
C             DMIN
C             ARBAYS
C             ARCOEF
C             BAYSWT
C             REDUCT
C             SDCOMP
C       ----------------------------------------------------------------
C
C     INPUTS:
C        SETX:  EXTERNAL SUBROUTINE DESIGNATION
C        Z:     ORIGINAL DATA, OUTPUT OF SUBROUTINE REDATA
C        X:     WORKING AREA
C        D:     WORKING AREA
C        LAG:   UPPER LIMIT OF ORDER OF AR-MODEL
C        N0:    INDEX OF THE END POINT OF THE FORMER SPAN
C        NS:    LENGTH OF BASIC LOCAL SPAN
C        KMAX:  MAXIMUM NUMBER OF PRECEDING MODELS STORED
C        KSW:   =0  CONSTANT VECTOR IS NOT INCLUDED AS A REGRESSOR
C               =1  CONSTANT VECTOR IS INCLUDED AS THE FIRST REGRESSOR
C        ISW:   PRINT OUT CONTROL
C        TITLE: TITLE OF DATA
C        MJ1:   ABSOLUTE DIMENSION OF X
C
C     OUTPUTS:
C        A:     AR-COEFFICIENTS OF THE CURRENT MODEL
C        SD:    INNOVATION VARIANCE OF THE CURRENT MODEL
C
cxx      IMPLICIT  REAL*8  ( A-H,O-Z )
CC      REAL*4  Z , TITLE
cc      REAL*4  TITLE
cc      DIMENSION  X(MJ1,1) , D(1) , A(1) , Z(1) , TITLE(1)
cc      DIMENSION  F(50,20) , AIC(21) , C(50) , B(50)
cx      DIMENSION  X(MJ1,1) , D(LAG+KSW+1), A(1) , Z(1)
cxx      DIMENSION  X(MJ1,1) , D(LAG+KSW+1), A(LAG+KSW) , Z(N)
cxx      DIMENSION  F(LAG+KSW,KMAX1) , AIC(KMAX1) , C(KMAX1) , B(LAG+KSW)
cxx      DIMENSION  SDD(LAG+KSW+1), AICC(LAG+KSW+1), DIC(LAG+KSW+1)
cxx      DIMENSION  B1(LAG+KSW), W(LAG+KSW+1)
      INTEGER :: N, LAG, N0, NS, KMAX1, KSW, ISW, MJ1, KC
      REAL(8) :: Z(N), X(MJ1,1), F(LAG+KSW,KMAX1), AIC(KMAX1),
     1           C(KMAX1), B(LAG+KSW), A(LAG+KSW), SD
      REAL(8) :: D(LAG+KSW+1), SDD(LAG+KSW+1), AICC(LAG+KSW+1),
     1           DIC(LAG+KSW+1), B1(LAG+KSW), W(LAG+KSW+1),
     2           AICB, AICM, SDMIN, PN
cc      DATA  KC / 0 /
      EXTERNAL  SETX
C
      K = LAG + KSW
      KMAX = KMAX1-1
C
C     ---  HOUSEHOLDER REDUCTION  ---
cc      CALL  REDUCT( SETX,Z,D,NS,N0,K,MJ1,LAG,X )
      CALL  REDUCT( SETX,Z,NS,N0,K,MJ1,LAG,X )
C
C
C     ---  BAYESIAN MODEL FITTED TO THE NEW SPAN  ---
cc      CALL  ARBAYS( X,D,K,LAG,NS,ISW,TITLE,MJ1,A,B,SD,AICB )
cxx      CALL ARBAYS( X,D,K,LAG,NS,ISW,MJ1,SDD,AICC,DIC,AICM,SDMIN,IMIN,
      CALL ARBAYS( X,D,K,NS,ISW,MJ1,SDD,AICC,DIC,AICM,SDMIN,IMIN,
     *             A,B1,B,W,SD,PN,AICB )
C
      IF( KC .EQ. 0 )  GO TO 110
C
C
C     ---  PREDICTION ERROR VARIANCE AND AIC OF THE FORMER MODELS  ---
      AIC(1) = AICB
      DO 30  J=1,KC
         DO 20  I=1,K
cxx   20    D(I) = F(I,J)
         D(I) = F(I,J)
   20    CONTINUE
        CALL  ARCOEF( D,K,A )
cc        CALL  SDCOMP( X,A,D,NS,K,MJ1,SD )
        CALL  SDCOMP( X,A,NS,K,MJ1,SD )
cxx   30 AIC(J+1) = NS*DLOG( SD ) + 2.D0
      AIC(J+1) = NS*DLOG( SD ) + 2.D0
   30 CONTINUE
C
C
C     ---  BAYESIAN WEIGHTS OF THE MODEL  ---
c-------------------------------   06/11/01
ccx      AICM = DMIN( AIC,KC )
      AICM = AIC(1)
      DO 33  I=1,KC
cxx   33 IF( AIC(I) .LT. AICM )  AICM = AIC(I)
      IF( AIC(I) .LT. AICM )  AICM = AIC(I)
   33 CONTINUE
c-------------------------------
      CALL  BAYSWT( AIC,AICM,KC,2,C )
C
      KC1 = KC+1
cc      WRITE( 6,3 )     C(1) , AIC(1)
cc      DO 35  I=2,KC1
cc         IM1 = I-1
cc   35 WRITE( 6,4 )     IM1 , C(I) , AIC(I)
C
C
C     ---  AVERAGING OF THE MODELS  ---
      DO 40  I=1,K
cxx   40 B(I) = B(I)*C(1)
      B(I) = B(I)*C(1)
   40 CONTINUE
      DO 70  J=1,KC
         DO 50  I=1,K
cxx   50    A(I) = F(I,J)
         A(I) = F(I,J)
   50    CONTINUE
         DO 60  I=1,K
cxx   60    B(I) = B(I) + C(J+1)*A(I)
         B(I) = B(I) + C(J+1)*A(I)
   60    CONTINUE
   70 CONTINUE
cc      WRITE( 6,5 )     (B(I),I=1,K)
C
C
C     ---  AR-COEFFICIENTS OF THE CURRENT BAYESIAN MODEL  ---
      CALL  ARCOEF( B,K,A )
C
C
C     ---  "PARCOR'S" STORED  ---
cxx      DO 100  J=1,KC
      DO 101  J=1,KC
         II = KC1-J
         DO 100  I=1,K
cxx  100 F(I,II+1) = F(I,II)
         F(I,II+1) = F(I,II)
  100    CONTINUE
  101 CONTINUE
  110 CONTINUE
      DO 120  I=1,K
cxx  120 F(I,1) = B(I)
      F(I,1) = B(I)
  120 CONTINUE
      KC = MIN0( KC+1,KMAX )
C
C
C     ---  PREDICTION ERROR VARIANCE COMPUTED  ---
cc      CALL  SDCOMP( X,A,D,NS,K,MJ1,SD )
      CALL  SDCOMP( X,A,NS,K,MJ1,SD )
C
      RETURN
C
cxx    3 FORMAT( ///1H ,13X,'AR-MODEL FITTED TO  !  BAYESIAN WEIGHTS  ! AIC
cxx     1 WITH RESPECT TO THE PRESENT DATA',/,10X,83(1H-),/,1H ,11X,'CURREN
cxx     2T BLOCK',9X,'!',F13.5,7X,'!',F21.3 )
cxx    4 FORMAT( 1H ,6X,I5,' PERIOD FORMER BLOCK  !',F13.5,7X,'!',F21.3 )
cxx    5 FORMAT( //1H ,'PARTIAL AUTOCORRELATION  B(I) (I=1,K)',/,(1X,10D13.
cxx     15) )
C
      END


text:
!> \file kronecker.f90
!! |brief Module defining Kronecker operations for matrix computation
!!
!! @author Innocent Souopgui
!!
!<

!> Module defining Kronecker product for matrix computation
!!
!<
module kronecker
  USE general_constant
  USE debug_tools

  private kron_product_check, kron_product_check_diag, kron_product_check_full
  !>Forms the Kronecker product of two matrices
  INTERFACE kron_product
    MODULE PROCEDURE kron_product_full
    MODULE PROCEDURE kron_product_diag
  END INTERFACE kron_product

  !>Check the input for the Kronecker product of two matrix
  INTERFACE kron_product_check
    MODULE PROCEDURE kron_product_check_full
    MODULE PROCEDURE kron_product_check_diag
  END INTERFACE kron_product_check

  !> matrix vector multiplication with the matrix given by a Kronecker product of 2/3 matrices
  INTERFACE kron_matvec
    MODULE PROCEDURE kron_matvec_full ! matrices are full and represented by 2d arrays
    MODULE PROCEDURE kron_matvec_diag ! matrices are diagonal and represented by vectors
    MODULE PROCEDURE kron_matvec_full3 ! matrices are full and represented by 2d arrays
    MODULE PROCEDURE kron_matvec_diag3 ! matrices are diagonal and represented by vectors
  END INTERFACE kron_matvec

CONTAINS

  !> \brief Computes the matrix vector product with the matrix given by a Kronecker product of two matrices.
  !! param[in] A first matrix of the kronecker product, size nrA x ncA
  !! param[in] B second matrix of the kronecker product, size nrB x ncB
  !! param[in] f input vector of size ncA*ncB
  !! param[out] g output vector of size nrA*nrB
  !! param[in] trans (optional) says if the transpose of the input matrices should be used instead of those matrices themselves.
  !! computes g = Kron(A,B)f if trans is not present or is set to TRUE, g = Kron(A,B)^Tf = (f^TKron(A,B) )^T if trans is present and set to TRUE
  !<
  SUBROUTINE kron_matvec_full(A, B, f, g, trans)
    REAL(cp), DIMENSION(:,:), INTENT(IN) :: A, B
    REAL(cp), DIMENSION(:), INTENT(IN)  :: f
    REAL(cp), DIMENSION(:), INTENT(OUT) :: g
    LOGICAL , INTENT(IN), OPTIONAL :: trans
    !local variables
    REAL(cp), DIMENSION(SIZE(B,1)) :: td!direct computation
    REAL(cp), DIMENSION(SIZE(B,2)) :: tt!transpose computation
    INTEGER :: nrA, ncA, nrB, ncB!number of row of X(nrX), !number of column of X(ncX)
    INTEGER :: iA, jA, if_start, if_end, ig_start, ig_end, nf_predicted, ng_predicted
    LOGICAL ll_trans

    IF( PRESENT(trans) )THEN
      ll_trans = trans
    ELSE
      ll_trans = .FALSE.
    END IF

    IF(ll_trans)THEN
      nf_predicted = SIZE(A,1)*SIZE(B,1)
      ng_predicted = SIZE(A,2)*SIZE(B,2)
    ELSE
      nf_predicted = SIZE(A,2)*SIZE(B,2)
      ng_predicted = SIZE(A,1)*SIZE(B,1)
    END IF
    IF( (SIZE(f)/=nf_predicted).OR.(SIZE(g)/=ng_predicted) )THEN
      CALL debug('Expecting data with following specification, based on the shape of A and B', tag=dALLWAYS)
      CALL debug(nf_predicted, 'f of size :', tag=dALLWAYS)
      CALL debug(ng_predicted, 'g of size :', tag=dALLWAYS)
      CALL debug('For debuging purpose, here are the shapes of Ux and Uy', tag=dALLWAYS)
      CALL debug(SHAPE(A), 'SHAPE(A)', tag=dALLWAYS)
      CALL debug(SHAPE(B), 'SHAPE(B)', tag=dALLWAYS)
      CALL stop_program('stopped from kron_matvec (with full matrices)')
    END IF

    nrA = SIZE(A,1)
    ncA = SIZE(A,2)
    nrB = SIZE(B,1)
    ncB = SIZE(B,2)

    g=0.0_cp
    IF(ll_trans)THEN! used the transposed
      if_start = 1
      if_end   = nrB
      DO iA = 1, nrA
        tt = MATMUL(f(if_start:if_end), B)
        ig_start = 1
        ig_end   = ncB
        DO jA=1,ncA
          g(ig_start:ig_end) = g(ig_start:ig_end) + A(iA,jA)*tt
          ig_start = ig_start + ncB
          ig_end   = ig_end   + ncB
        END DO
        if_start = if_start + nrB
        if_end   = if_end   + nrB
      END DO
    ELSE
      if_start = 1
      if_end   = ncB
      DO jA = 1, ncA
        td = MATMUL(B, f(if_start:if_end))
        ig_start = 1
        ig_end   = nrB
        DO iA=1,nrA
          g(ig_start:ig_end) = g(ig_start:ig_end) + A(iA,jA)*td
          ig_start = ig_start + nrB
          ig_end   = ig_end   + nrB
        END DO
        if_start = if_start + ncB
        if_end   = if_end   + ncB
      END DO
    END IF
  END SUBROUTINE kron_matvec_full

  !> @brief Computes the matrix vector product with the matrix given
  !!   by a Kronecker product of 3 matrices.
  !! @param[in] A first matrix of the kronecker product, size nrA x ncA
  !! @param[in] B second matrix of the kronecker product, size nrB x ncB
  !! @param[in] C third matrix of the kronecker product, size nrC x ncC
  !! @param[in] f input vector of size ncA*ncB*ncC
  !! @param[out] g output vector of size nrA*nrB*nrC
  !! @param[in] trans (optional) says if the transpose of the input matrices
  !!   should be used instead of those matrices themselves.
  !! @details computes g = Kron(A,B,C)f if trans is not present or is set to TRUE,
  !!   g = Kron(A,B)^Tf = (f^TKron(A,B,C) )^T if trans is present and set to TRUE
  !<
  subroutine kron_matvec_full3(A, B, C, f, g, trans)
    real(cp), dimension(:,:), intent(in) :: A, B, C
    real(cp), dimension(:), intent(in)  :: f
    real(cp), dimension(:), intent(out) :: g
    logical , intent(in), optional :: trans
    !local variables
    real(cp), dimension(size(B,1)*size(C,1)) :: td!direct computation
    real(cp), dimension(size(B,2)*size(C,2)) :: tt!transpose computation
    integer :: nrA, ncA, nrB, ncB, nrC, ncC!number of row of X(nrX), !number of column of X(ncX)
    integer :: iA, jA, iB, jB, if_start, if_end, ig_start&
             , ig_end, nf_predicted, ng_predicted
    logical ll_trans

    if( present(trans) )then
      ll_trans = trans
    else
      ll_trans = .false.
    end if

    if(ll_trans)then
      nf_predicted = size(A,1)*size(B,1)*size(C,1)
      ng_predicted = size(A,2)*size(B,2)*size(C,2)
    else
      nf_predicted = size(A,2)*size(B,2)*size(C,2)
      ng_predicted = size(A,1)*size(B,1)*size(C,1)
    end if
    if( (size(f)/=nf_predicted).or.(size(g)/=ng_predicted) )then
      call debug('Expecting data with following specification, based on the shape of A, B, C', tag=dALLWAYS)
      call debug(nf_predicted, 'f of size :', tag=dALLWAYS)
      call debug(ng_predicted, 'g of size :', tag=dALLWAYS)
      call debug('For debuging purpose, here are the shapes of A, B, C', tag=dALLWAYS)
      call debug(shape(A), 'SHAPE(A)', tag=dALLWAYS)
      call debug(shape(B), 'SHAPE(B)', tag=dALLWAYS)
      call debug(shape(C), 'SHAPE(C)', tag=dALLWAYS)
      call stop_program('stopped from kron_matvec (with full matrices)')
    end if

    nrA = size(A,1)
    ncA = size(A,2)
    nrB = size(B,1)
    ncB = size(B,2)
    nrC = size(C,1)
    ncC = size(C,2)

    g=0.0_cp
    if(ll_trans)then! used the transposed
      if_start = 1
      if_end   = nrB*nrC
      do iA = 1, nrA
        call kron_matvec(B, C, f(if_start:if_end), tt, trans=.true.)
        ig_start = 1
        ig_end   = ncB*ncC
        do jA=1,ncA
          g(ig_start:ig_end) = g(ig_start:ig_end) + A(iA,jA)*tt
          ig_start = ig_start + ncB*ncC
          ig_end   = ig_end   + ncB*ncC
        end do
        if_start = if_start + nrB*nrC
        if_end   = if_end   + nrB*nrC
      end do
    else
      if_start = 1
      if_end   = ncB*ncC
      do jA = 1, ncA
        call kron_matvec(B, C, f(if_start:if_end), td, trans=.false.)
        ig_start = 1
        ig_end   = nrB*nrC
        do iA=1,nrA
          g(ig_start:ig_end) = g(ig_start:ig_end) + A(iA,jA)*td
          ig_start = ig_start + nrB*nrC
          ig_end   = ig_end   + nrB*nrC
        end do
        if_start = if_start + ncB*ncC
        if_end   = if_end   + ncB*ncC
      end do
    end if
  end subroutine kron_matvec_full3

  !> \brief Computes the matrix vector product with the matrix given by
  !!   a Kronecker product of two matrices
  !! the difference with kron_matvec_full is that the matrices are diagonal
  !!   and represented by vectors of diagonal elements.
  !! @see kron_matvec_full for the description of parameters
  !<
  SUBROUTINE kron_matvec_diag(A, B, f, g)
    REAL(cp), DIMENSION(:), INTENT(IN) :: A, B, f
    REAL(cp), DIMENSION(:), INTENT(OUT) :: g
    !local variables
    INTEGER :: nA, nB !number of elements of X
    INTEGER :: iA, i_start, i_end, n_predicted

    nA = SIZE(A)
    nB = SIZE(B)
    n_predicted = nA*nB
    IF( (SIZE(f)/=n_predicted).OR.(SIZE(g)/=n_predicted) )THEN
      CALL debug('Expecting data with following specification, based on the shape of A and B', tag=dALLWAYS)
      CALL debug(n_predicted, 'f and g of size :', tag=dALLWAYS)
      CALL debug('For debuging purpose, here are the SIZE of A and B', tag=dALLWAYS)
      CALL debug(SIZE(A), 'SIZE(A)', tag=dALLWAYS)
      CALL debug(SIZE(B), 'SIZE(B)', tag=dALLWAYS)
      CALL stop_program('stopped from kron_matvec_diag')
    END IF

    g=0.0_cp

    i_start = 1
    i_end   = nB
    DO iA = 1, nA
      g(i_start:i_end) = A(iA)*( B*f(i_start:i_end) )
      i_start = i_start + nB
      i_end   = i_end   + nB
    END DO
  END SUBROUTINE kron_matvec_diag

  !> \brief Computes the matrix vector product with the matrix given by
  !!   a Kronecker product of 3 matrices
  !! the difference with kron_matvec_full3 is that the matrices are diagonal
  !!   and represented by vectors of diagonal elements.
  !! @see kron_matvec_full3 for the description of parameters
  !<
  subroutine kron_matvec_diag3(A, B, C, f, g)
    real(cp), dimension(:), intent(in) :: A, B, C, f
    real(cp), dimension(:), intent(out) :: g
    !local variables
    integer :: nA, nB, nC !number of elements of X
    integer :: iA, iB, i_start, i_end, n_predicted

    nA = size(A)
    nB = size(B)
    nC = size(C)
    n_predicted = nA*nB*nC
    if( (size(f)/=n_predicted).or.(size(g)/=n_predicted) )then
      call debug('Expecting data with following specification, based on the shape of A, B, C', tag=dALLWAYS)
      call debug(n_predicted, 'f and g of size :', tag=dALLWAYS)
      call debug('For debuging purpose, here are the SIZE of A and B', tag=dALLWAYS)
      call debug(size(A), 'size(A)', tag=dALLWAYS)
      call debug(size(B), 'size(B)', tag=dALLWAYS)
      call debug(size(C), 'size(C)', tag=dALLWAYS)
      call stop_program('stopped from kron_matvec_diag3')
    end if

    g=0.0_cp

    i_start = 1
    i_end   = nC
    do iA = 1, nA
        do iB = 1, nB
            g(i_start:i_end) = A(iA)*B(iB)*( C*f(i_start:i_end) )
            i_start = i_start + nC
            i_end   = i_end   + nC
        end do
    end do
  end subroutine kron_matvec_diag3

  !> Checks parameters of the kronecker product
  !! @see kron_product_full
  SUBROUTINE kron_product_check_full(A, B, C)
    REAL(cp), DIMENSION(:,:), INTENT(IN) :: A, B, C
    !local variables
    INTEGER :: nrA, ncA, nrB, ncB!number of row of X(nrX), !number of column of X(ncX)

    nrA = SIZE(A,1)
    ncA = SIZE(A,2)
    nrB = SIZE(B,1)
    ncB = SIZE(B,2)
    IF( ( SIZE(C,1)/=nrA*nrB ).OR.( ( SIZE(C,2)/=ncA*ncB ) ) )THEN
      CALL debug('Expecting C with following specification, based on the shape of A and B', tag=dALLWAYS)
      CALL debug((/nrA*nrB, ncA*ncB/), 'expected SHAPE(C)', tag=dALLWAYS)
      CALL debug('Get the following, check and try aigain', tag=dALLWAYS)
      CALL debug(SHAPE(C), 'SHAPE(C)', tag=dALLWAYS)
      CALL debug('For debuging purpose, here are the shapes of A and B', tag=dALLWAYS)
      CALL debug(SHAPE(A), 'SHAPE(A)', tag=dALLWAYS)
      CALL debug(SHAPE(B), 'SHAPE(B)', tag=dALLWAYS)
      CALL stop_program('Bad shape for output matrix in kronecker product')
    END IF
  END SUBROUTINE kron_product_check_full

  !> Checks parameters of the kronecker product
  !! @see kron_product_full
  SUBROUTINE kron_product_check_diag(A, B, C)
    REAL(cp), DIMENSION(:), INTENT(IN) :: A, B, C
    !local variables
    INTEGER :: nA, nB

    nA = SIZE(A)
    nB = SIZE(B)
    IF( SIZE(C)/=nA*nB )THEN
      CALL debug('Expecting C with following specification, based on the SIZE of A and B', tag=dALLWAYS)
      CALL debug(nA*nB, 'expected SIZE(C)', tag=dALLWAYS)
      CALL debug('Get the following, check and try aigain', tag=dALLWAYS)
      CALL debug(SIZE(C), 'SIZE(C)', tag=dALLWAYS)
      CALL debug('For debuging purpose, here are the SIZEs of A and B', tag=dALLWAYS)
      CALL debug(SIZE(A), 'SIZE(A)', tag=dALLWAYS)
      CALL debug(SIZE(B), 'SIZE(B)', tag=dALLWAYS)
      CALL stop_program('Bad SIZE for output matrix in diagonal kronecker product')
    END IF
  END SUBROUTINE kron_product_check_diag

  !> \brief Forms the Kronecker product of two matrices
  !! \param[in] A first matrix in the Kronecker product
  !! \param[in] B second matrix in the Kronecker product
  !! \param[out] C resulting matrix
  !<
  SUBROUTINE kron_product_full(A, B, C)
    REAL(cp), DIMENSION(:,:), INTENT(IN) :: A, B
    REAL(cp), DIMENSION(:,:), INTENT(OUT) :: C
    !local variables
    INTEGER :: nrA, ncA, nrB, ncB, nrC, ncC!number of row of X(nrX), !number of column of X(ncX)
    INTEGER :: iA, jA, jB, iC_start, iC_end, jC

    !checking the input parameters
    CALL kron_product_check(A, B, C)
    nrA = SIZE(A,1)
    ncA = SIZE(A,2)
    nrB = SIZE(B,1)
    ncB = SIZE(B,2)
    nrC = SIZE(C,1)
    ncC = SIZE(C,2)

    DO jA=1,ncA
      !computing a bloc of columns of C
      jC = (jA-1)*ncB + 1
      DO jB=1,ncB
        !computing a column of C
        iC_start = 1
        iC_end   = nrB
        DO iA=1,nrA
          C(ic_start:ic_end, jC) = A(iA,jA)*B(:,jB)
          iC_start = iC_start + nrB
          iC_end   = iC_end   + nrB
        END DO
        jC = jC + 1
      END DO
    END DO
  END SUBROUTINE kron_product_full

  !> \brief Forms the Kronecker product of two matrices
  !! the difference with kron_matvec_full is that the matrices are diagonal and represented by vectors of diagonal elements.
  !! @see kron_product_full for the description of parameters
  !<
  SUBROUTINE kron_product_diag(A, B, C)
    REAL(cp), DIMENSION(:), INTENT(IN) :: A, B
    REAL(cp), DIMENSION(:), INTENT(OUT) :: C
    !local variables
    INTEGER :: nA, nB
    INTEGER :: iA, iC_start, iC_end

    !checking the input parameters
    CALL kron_product_check(A, B, C)
    nA = SIZE(A)
    nB = SIZE(B)

    iC_start = 1
    iC_end   = nB
    DO iA=1,nA
      C(ic_start:ic_end) = A(iA)*B
      iC_start = iC_start + nB
      iC_end   = iC_end   + nB
    END DO
  END SUBROUTINE kron_product_diag

end module kronecker

max_stars_repo_name: souopgui/da-toolbox

max_issues_repo_name: souopgui/da-toolbox

max_forks_repo_name: souopgui/da-toolbox