text:
      SUBROUTINE SBDT03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
     $                   RESID )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            KD, LDU, LDVT, N
      REAL               RESID
*     ..
*     .. Array Arguments ..
      REAL               D( * ), E( * ), S( * ), U( LDU, * ),
     $                   VT( LDVT, * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  SBDT03 reconstructs a bidiagonal matrix B from its SVD:
*     S = U' * B * V
*  where U and V are orthogonal matrices and S is diagonal.
*
*  The test ratio to test the singular value decomposition is
*     RESID = norm( B - U * S * VT ) / ( n * norm(B) * EPS )
*  where VT = V' and EPS is the machine precision.
*
*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          Specifies whether the matrix B is upper or lower bidiagonal.
*          = 'U':  Upper bidiagonal
*          = 'L':  Lower bidiagonal
*
*  N       (input) INTEGER
*          The order of the matrix B.
*
*  KD      (input) INTEGER
*          The bandwidth of the bidiagonal matrix B.  If KD = 1, the
*          matrix B is bidiagonal, and if KD = 0, B is diagonal and E is
*          not referenced.  If KD is greater than 1, it is assumed to be
*          1, and if KD is less than 0, it is assumed to be 0.
*
*  D       (input) REAL array, dimension (N)
*          The n diagonal elements of the bidiagonal matrix B.
*
*  E       (input) REAL array, dimension (N-1)
*          The (n-1) superdiagonal elements of the bidiagonal matrix B
*          if UPLO = 'U', or the (n-1) subdiagonal elements of B if
*          UPLO = 'L'.
*
*  U       (input) REAL array, dimension (LDU,N)
*          The n by n orthogonal matrix U in the reduction B = U'*A*P.
*
*  LDU     (input) INTEGER
*          The leading dimension of the array U.  LDU >= max(1,N)
*
*  S       (input) REAL array, dimension (N)
*          The singular values from the SVD of B, sorted in decreasing
*          order.
*
*  VT      (input) REAL array, dimension (LDVT,N)
*          The n by n orthogonal matrix V' in the reduction
*          B = U * S * V'.
*
*  LDVT    (input) INTEGER
*          The leading dimension of the array VT.
*
*  WORK    (workspace) REAL array, dimension (2*N)
*
*  RESID   (output) REAL
*          The test ratio:  norm(B - U * S * V') / ( n * norm(A) * EPS )
*
* ======================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
      REAL               BNORM, EPS
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ISAMAX
      REAL               SASUM, SLAMCH
      EXTERNAL           LSAME, ISAMAX, SASUM, SLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEMV
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN, REAL
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      RESID = ZERO
      IF( N.LE.0 )
     $   RETURN
*
*     Compute B - U * S * V' one column at a time.
*
      BNORM = ZERO
      IF( KD.GE.1 ) THEN
*
*        B is bidiagonal.
*
         IF( LSAME( UPLO, 'U' ) ) THEN
*
*           B is upper bidiagonal.
*
            DO 20 J = 1, N
               DO 10 I = 1, N
                  WORK( N+I ) = S( I )*VT( I, J )
   10          CONTINUE
               CALL SGEMV( 'No transpose', N, N, -ONE, U, LDU,
     $                     WORK( N+1 ), 1, ZERO, WORK, 1 )
               WORK( J ) = WORK( J ) + D( J )
               IF( J.GT.1 ) THEN
                  WORK( J-1 ) = WORK( J-1 ) + E( J-1 )
                  BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J-1 ) ) )
               ELSE
                  BNORM = MAX( BNORM, ABS( D( J ) ) )
               END IF
               RESID = MAX( RESID, SASUM( N, WORK, 1 ) )
   20       CONTINUE
         ELSE
*
*           B is lower bidiagonal.
*
            DO 40 J = 1, N
               DO 30 I = 1, N
                  WORK( N+I ) = S( I )*VT( I, J )
   30          CONTINUE
               CALL SGEMV( 'No transpose', N, N, -ONE, U, LDU,
     $                     WORK( N+1 ), 1, ZERO, WORK, 1 )
               WORK( J ) = WORK( J ) + D( J )
               IF( J.LT.N ) THEN
                  WORK( J+1 ) = WORK( J+1 ) + E( J )
                  BNORM = MAX( BNORM, ABS( D( J ) )+ABS( E( J ) ) )
               ELSE
                  BNORM = MAX( BNORM, ABS( D( J ) ) )
               END IF
               RESID = MAX( RESID, SASUM( N, WORK, 1 ) )
   40       CONTINUE
         END IF
      ELSE
*
*        B is diagonal.
*
         DO 60 J = 1, N
            DO 50 I = 1, N
               WORK( N+I ) = S( I )*VT( I, J )
   50       CONTINUE
            CALL SGEMV( 'No transpose', N, N, -ONE, U, LDU, WORK( N+1 ),
     $                  1, ZERO, WORK, 1 )
            WORK( J ) = WORK( J ) + D( J )
            RESID = MAX( RESID, SASUM( N, WORK, 1 ) )
   60    CONTINUE
         J = ISAMAX( N, D, 1 )
         BNORM = ABS( D( J ) )
      END IF
*
*     Compute norm(B - U * S * V') / ( n * norm(B) * EPS )
*
      EPS = SLAMCH( 'Precision' )
*
      IF( BNORM.LE.ZERO ) THEN
         IF( RESID.NE.ZERO )
     $      RESID = ONE / EPS
      ELSE
         IF( BNORM.GE.RESID ) THEN
            RESID = ( RESID / BNORM ) / ( REAL( N )*EPS )
         ELSE
            IF( BNORM.LT.ONE ) THEN
               RESID = ( MIN( RESID, REAL( N )*BNORM ) / BNORM ) /
     $                 ( REAL( N )*EPS )
            ELSE
               RESID = MIN( RESID / BNORM, REAL( N ) ) /
     $                 ( REAL( N )*EPS )
            END IF
         END IF
      END IF
*
      RETURN
*
*     End of SBDT03
*
      END

max_stars_repo_name: nashmit/lapack-release

max_issues_repo_name: nashmit/lapack-release

max_forks_repo_name: nashmit/lapack-release


text:
c***********************************************************************
c  READP_F: main entry point for Read_Phoenics module
c  Author: David Knight, AVS(UK) Ltd
c  Creation Date: 26.8.97
c  Notes:
c***********************************************************************
      integer function READP_F(OBJ_ID, MASK, SEQ_NUM)
      include 'avs/omf.inc'
C **
C PARAMETERS DEFINING INTERNAL ARRAY SIZES, THESE MAY NEED TO BE
C TO BE INCREASED FOR LARGE PROBLEMS. HERE 'BIG GRID' MEANS THE
C ENTIRE MULTIBLOCK MESH
C **
c .. Maximum number of bodies
      PARAMETER(IBMAX=100)
c .. Maximum number of links
      PARAMETER(NLMAX=100)
c .. Maximum number of nodes in the big grid
      PARAMETER(NODMAX=100000)
c .. Maximum number of elements in the big grid
      PARAMETER(IELMAX=100000)
c .. Maximum number of nodal results (multiple of NODMAX)
      PARAMETER(NRMAX = NODMAX * 2)
c .. If NRMAX*3 + IELMAX is greater than 8*IELMAX IWORK must be
c .. dimensioned to   NRMAX*3 + IELMAX
      LOGICAL XDIR,PDIR,BFC,CARTES,LOG(50),NOBLK,COMM,CCC1
      LOGICAL ONEPHS,CRT1,CCC2,CRT2,FVPOR,FPRPS
      CHARACTER*4 NAME(50)
      CHARACTER*40 TITLE
      CHARACTER*80 PHINAM,XYZNAM,PATNAM
      CHARACTER*40 MESSAGE
      INTEGER IJKBLK(6,IBMAX),NMAP(NODMAX),IWORK(8*IELMAX)
      INTEGER LINK(11,NLMAX),LIMBLK(3,IBMAX),NRESV(50)
      REAL WORK(3*NODMAX + 3*NRMAX),XFRAC(256),TEMP(NODMAX)
c .. express objects
      INTEGER OBJ_ID(OIDSIZ), FLD_ID(OIDSIZ), CDA_ID(OIDSIZ),
     .	NDA_ID(OIDSIZ)
      INTEGER MASK, SEQ_NUM
c .. IO units
      IOPAT = 50
c .. end NEP stuff to stdout
      IONEP = 51
      IOPHI = 52
      IOPHDA = 53
      IOGRD = 54
      IOPRE = 55
      IOCON = 56
      COMM = .TRUE.
      MESSAGE = ' '
c .. get the input parameters
      CALL GET_PARAMS(OBJ_ID,PHINAM,XYZNAM,PATNAM,
     .	PDIR,XDIR,FLD_ID)
c .. assumed size of reals and ints (increase for 64bit)
      LBYTE = 4
      LENREF = 256
c .. open the NEP file
      CALL ONEPNP(IONEP,IERR)
      IF (IERR.GT.0) GOTO 95
c .. open the PATGEO file, if there isn't one, then use
c .. default parameters from the PHI file
      IF(COMM) WRITE(*,*)'Opening PATGEO file'
      CALL OPATNP(IONEP,IOPAT,PATNAM,IERR)
c .. no patgeo file
      IF (IERR.GT.0) THEN
        WRITE(*,*) 'PATGEO file not found. Using defaults.'
      ENDIF
      IERR = 0
c .. Open the PHI file, create a temp DA file if required
      WRITE(*,*)'Opening PHI file'
      IF(PDIR) THEN
        LEN = LBYTE*LENREF
c ..   Open with correct record length
         OPEN(IOPHDA,FILE=PHINAM,STATUS='OLD',RECL=LEN,
     +        ACCESS='DIRECT',ERR=97)
c ..   ASCII file
      ELSE
c ..   Open the PHI file
        OPEN(IOPHI,FILE=PHINAM,STATUS='OLD',ERR=97)
        LEN = LBYTE*LENREF
        OPEN(IOPHDA,STATUS='SCRATCH',RECL=LEN,
     +      ACCESS='DIRECT',ERR=97)
c ..   Convert to DA
        CALL CPHINP(IONEP,IOPHI,IOPHDA,IWORK,WORK,LENREF,IERR)
        IF(IERR .GT. 0) GOTO 97
      ENDIF
c .. read and output the file title
      CALL RPHINP(IOPHDA,TITLE,NX,NY,NZ,LOG,NAME,CARTES,ONEPHS,BFC,
     .            IONEP,IERR)
      IF(IERR .GT. 0) GOTO 99
      WRITE(IONEP,'(A40)') TITLE
c .. Set the title string
      CALL WRITE_TITLE(OBJ_ID,TITLE)
c .. size of the mesh
      NEBIG = NX*NY*NZ
      NNBIG = (NX+1)*(NY+1)*(NZ+1)
c .. read the block and inter block patches
      if (COMM) WRITE(*,*)'Reading PATGEO file'
      call RPATNP(IONEP,IOPAT,IBMAX,NODMAX,NLMAX,NX,NY,NZ,IJKBLK,NMAP,
     +             NBLOCK,NUMNOD,LINK,NLINK,IERR)
      IF(IERR .GT. 0) GOTO 95
c .. write the block definitions
c .. get the model coordinates
      if (COMM) write(*,*) 'Reading the model coordinates'
      if (BFC) then
      	if (XDIR) then
          LEN = LBYTE*LENREF
          OPEN(UNIT=IOGRD,FILE=XYZNAM,STATUS='OLD',ERR=92,RECL=LEN,
     +                   ACCESS='DIRECT')
          LENX = LENREF
        ELSE
          OPEN(UNIT=IOGRD,FILE=XYZNAM,STATUS='OLD',ERR=92)
          LENX = 5
        ENDIF
        CALL RNODNP(IONEP,IOGRD,XDIR,IJKBLK,NBLOCK,LENX,
     +           NODMAX,NX,NY,NZ,WORK(3*NODMAX+1),WORK,IERR)
        IF(IERR .GT. 0) GOTO 92
      ELSE
        KX = NX + 1
        KY = NY + 1
        KZ = NZ + 1
        IWOFF = 3*KX*KY*KZ + 1
        CALL CCALNP(IOPHDA,CARTES,NX,NY,NZ,XFRAC,WORK(IWOFF+LENREF),
     +            WORK(2*LENREF+IWOFF),LENREF,NUMNOD,WORK,IONEP,IERR)
        IF(IERR .GT. 0) GOTO 92
      ENDIF
      IF(COMM) WRITE(*,*)'Writing the coordinates'
c .. write the nodes
c      CALL WNODNP(IONEP,NMAP,WORK,NUMNOD)
c .. write field coordinates
      CALL WRITE_COORDS(fld_id,NMAP,WORK,NUMNOD)
c .. Calculate the elements definitions
      IF(COMM) WRITE(*,*)'Calculating the element definitions'
      CALL CELENP(IONEP,NMAP,IJKBLK,NBLOCK,IELMAX,IWORK,
     +             NUMELE,IERR)
      IF(IERR .GT. 0) GOTO 91
c .. write the element definitions
      IF(COMM) WRITE(*,*)'Writing the element definitions'
c      CALL WELENP(IONEP,IWORK,NUMELE)
c .. write the element definitions into the field
      CALL WRITE_ELEMENTS(FLD_ID,IWORK,NUMELE)
c .. find the prps numbers
      IF(COMM) WRITE(*,*)'Reading the PRPS numbers'
      CALL PRPSNP(IONEP,IOPHDA,LENREF,LOG,NAME,IJKBLK,
     +     NBLOCK,NX,NY,NZ,WORK(2*NEBIG+1),NUMRES,IWORK,WORK,IERR)
      IF(IERR .GT. 0) GOTO 91
c .. calculate the number of results to output
c .. here, the nnode_data = numres and ncell_data = numres + 1
c .. unless the PRPS or VPOR are already included
      NRES = NUMRES * 2
      FPRPS = .FALSE.
      FVPOR = .FALSE.
      DO 15 IR = 1,NUMRES
        IF(NAME(IR).EQ.'PRPS' )FPRPS = .TRUE.
        IF(NAME(IR).EQ.'VPOR') FVPOR = .TRUE.
 15   CONTINUE
      IF(FPRPS .AND. FVPOR) THEN
        NRES = NRES - 2
      ELSEIF(FPRPS) THEN
        NRES = NRES - 1
      ELSEIF(FVPOR) THEN
        NRES = NRES
      ELSE
        NRES = NRES + 1
      ENDIF
      WRITE(IONEP,*)NRES
c .. initialise field data
      CALL INIT_FLDDATA(FLD_ID,NUMRES,CDA_ID,NDA_ID)
c .. write the prps numbers out
      IF(COMM) print*,'Writing the PRPS numbers'
      IEOFF = 0
      WRITE(IONEP,*)'PRPS ELEMC'
      DO 16 IB = 1,NBLOCK
        NE = (IJKBLK(2,IB)-IJKBLK(1,IB)+1)*
     +     (IJKBLK(4,IB)-IJKBLK(3,IB)+1)*(IJKBLK(6,IB)-IJKBLK(5,IB)+1)
c       WRITE(IONEP,*)NE
c       WRITE(IONEP,'(3(I8,1X,E13.6,1X))')
c     +      (IEOFF+IE,WORK(IEOFF+IE),IE=1,NE)
        IEOFF = IEOFF + NE
c .. write into the cell data array
      call WRITE_DATA(CDA_ID,0,'PRPS',NUMELE,WORK)
 16   CONTINUE
c .. find the velocity locations
      CALL FNDVNP(1,NAME,NUMRES,CCC1,CRT1,LOCU1,LOCV1,LOCW1)
      CALL FNDVNP(2,NAME,NUMRES,CCC2,CRT2,LOCU2,LOCV2,LOCW2)
c .. for each result stored, write into the cell and node data arrays
c .. we ignore the velocity components at this stage
      DO 100 IR = 1,NUMRES
C
        IF(NAME(IR) .NE. 'PRPS' .AND. IR.NE.LOCU1 .AND. IR.NE.LOCV1
     +      .AND. IR.NE.LOCW1 .AND. IR.NE.LOCU2 .AND. IR.NE.LOCV2
     +      .AND. IR.NE.LOCW2) THEN
          NOBLK = (NAME(IR).EQ.'TEM1' .OR. NAME(IR).EQ.'DEN1')
          IF(COMM) WRITE(*,*)'Reading ',NAME(IR)
          WRITE(IONEP,*)NAME(IR),' ELEMC'
c ..     Read the results
          CALL RDVLNP(IONEP,IOPHDA,IR,NUMRES,LENREF,IJKBLK,NBLOCK,
     +                NX,NY,NZ,WORK(NEBIG+1),WORK,IERR)
          IF(IERR .GT. 0) GOTO 91
c ..     Write the scalar results
          IF(COMM) WRITE(*,*)'Writing Cell centre ',NAME(IR)
          IEP = 0
          DO 17 IB = 1,NBLOCK
            CALL CSCANP(IB,IJKBLK,WORK,IWORK,NOBLK,NX,NY,NZ,
     +        WORK(NEBIG+1),IWORK(NEBIG+1),NRES)
c            WRITE(IONEP,*)NRES
c            WRITE(IONEP,'(3(I8,1X,E13.6,1X))')
c     +      (IWORK(NEBIG+I),WORK(NEBIG+I),I=1,NRES)
c .. write into an array
            DO 34 I = 1,NRES
              NE = IWORK(NEBIG+I)
      		TEMP(NE) = WORK(NEBIG+I)
 34         CONTINUE
 17       CONTINUE
c .. write into express object
      	CALL WRITE_DATA(CDA_ID,IR,NAME(IR),NUMELE,TEMP)
      	ICELL=ICELL+1
          IF( NAME(IR) .EQ. 'VPOR') GOTO 25
c ..
c ..     Nodally average the results
          IF(COMM) WRITE(*,*)'Writing Nodally Averaged ',NAME(IR)
          WRITE(IONEP,*)NAME(IR),' NODAL'
          IEOFF = 0
          IN = 0
          DO 20 IB = 1,NBLOCK
            CALL NAVGNP(IWORK,WORK,IJKBLK,IB,NOBLK,NMAP,NBLOCK,
     +        LINK,NLINK,NRMAX,WORK(NEBIG+1),IWORK(NEBIG+1),
     +         IWORK(NEBIG+NRMAX+1),NR,IONEP,IERR)
            IF(IERR .GT. 0) GOTO 91
c ..       Write the nodally averaged results
c            WRITE(IONEP,*)NR
c            WRITE(IONEP,'(3(I8,1X,I3,1X,E13.6,1X))')
c     +      (IWORK(NEBIG+NRMAX+I),IWORK(NEBIG+I),WORK(NEBIG+I),I=1,NR)
c .. write array
            DO 35 I = 1,NR
              NN = IWORK(NEBIG+NRMAX+I)
              TEMP(NN) = WORK(NEBIG+I)
 35	      CONTINUE
 20       CONTINUE
      	CALL WRITE_DATA(NDA_ID,IR-1,NAME(IR),NUMNOD,TEMP)
 25       CONTINUE
        ENDIF
 100  CONTINUE
c .. Write out velocities
      IF(COMM) WRITE(*,*)'Writing First Phase Velocities'
      CALL VELONP(LOCU1,LOCV1,LOCW1,IONEP,IOPHDA,LENREF,IJKBLK,NX,
     +    NY,NZ,NEBIG,1,NUMRES,IWORK,WORK,CCC1,CRT1,NBLOCK,NRMAX,NRESV,
     +    NUMNOD,NAME,NMAP,LINK,NLINK,IERR,
     +    NUMELE,TEMP,CDA_ID,NDA_ID,OIDSIZ)
      IF(IERR .GT. 0) GOTO 93
c .. Second phase velocity if required
      IF(.NOT. ONEPHS) THEN
        IF(COMM) WRITE(*,*)'Writing Second Phase Velocities'
        CALL VELONP(LOCU2,LOCV2,LOCW2,IONEP,IOPHDA,LENREF,IJKBLK,NX,
     +    NY,NZ,NEBIG,2,NUMRES,IWORK,WORK,CCC2,CRT2,NBLOCK,NMRES,NRESV,
     +    NUMNOD,NAME,NMAP,LINK,NLINK,IERR,
     +	NUMELE,TEMP,CDA_ID,NDA_ID,OIDSIZ)
        IF(IERR .GT. 0) GOTO 93
      ENDIF
c .. write links (unused at present)
c      IF(NLINK .GT. 0)CALL WRTLNK(IONEP,LINK,NLINK,IJKBLK,NBLOCK)
c .. normal exit
      readp_f = 1
      goto 90
c .. abnormal exits
98	MESSAGE = 'Error opening PATGEO file.'
      GOTO 99
97	MESSAGE = 'Error opening PHI file'
      GOTO 99
96	MESSAGE = 'Error opening XYZ file'
      GOTO 99
95	MESSAGE = 'Error opening NEP file'
      GOTO 99
94	MESSAGE = 'Error reading PATGEO file'
      GOTO 99
93	MESSAGE = 'Error writing velocities'
      GOTO 99
92	MESSAGE = 'Error getting coordinates'
      GOTO 99
91	MESSAGE = 'Error writing nodally averaged data'
99	readp_f = 0
90	CONTINUE
      WRITE(*,*) MESSAGE
      CALL WRITE_MESSAGE(OBJ_ID,MESSAGE)
c .. close all files
      CALL FCLSNP(IOPHI,IOPHDA,IONEP,IOPAT,IOGRD)
      RETURN
      END
      SUBROUTINE GET_PARAMS(OBJ_ID,PHINAM,XYZNAM,PATNAM,
     .    PDIR,XDIR,FLD_ID)
c**********************************************************************
c Get the input parameters from the Express object
c**********************************************************************
      INCLUDE 'avs/omf.inc'
      INTEGER OBJ_ID(OIDSIZ),FLD_ID(OIDSIZ),TMP_ID(OIDSIZ)
      CHARACTER*80 PHINAM,XYZNAM,PATNAM
      LOGICAL PDIR,XDIR
      INTEGER PDIR_ID(OIDSIZ), XDIR_ID(OIDSIZ)
      INTEGER ISTAT, PDIR_VAL, XDIR_VAL
c .. filename length
      NLEN = 80
c .. PHI filename (required)
c .. NB MS Fortran is really dopey, if a filename is blank
c .. it prompts the user for one, which isnt very helpful if
c .. its a Windows app. Hence the 'NULL' below
      ISTAT =  OMFget_name_str_val(OBJ_ID,
     .	OMFstr_to_name('phi_name'),
     .	PHINAM, NLEN)
      IF (ISTAT.LT.1.OR.PHINAM(1:2).EQ.' ') THEN
        PHINAM = 'NULL'
      ENDIF
c .. XYZ filename (optional)
      ISTAT =  OMFget_name_str_val(obj_id,
     .	OMFstr_to_name('xyz_name'),
     .	XYZNAM, NLEN)
      IF (ISTAT.LT.1.OR.XYZNAM(1:2).EQ.' ') THEN
        XYZNAM = 'NULL'
      ENDIF
c .. PATGEO filename (optional)
      ISTAT =  OMFget_name_str_val(obj_id,
     .	OMFstr_to_name('pat_name'),
     .	PATNAM, NLEN)
      IF (ISTAT.LT.1.OR.PATNAM(1:2).EQ.' ') THEN
        PATNAM = 'NULL'
      ENDIF
c .. phi binary ?
      ISTAT = OMFfind_subobj(obj_id, OMFstr_to_name('phi_bin'),
     .	OM_OBJ_RD, pdir_id)
      ISTAT = OMFget_int_val(pdir_id,pdir_val)
c .. xyz binary ?
      ISTAT = OMFfind_subobj(obj_id, OMFstr_to_name('xyz_bin'),
     .	OM_OBJ_RD, xdir_id)
      ISTAT = OMFget_int_val(xdir_id,xdir_val)
c .. set logicals
      pdir = pdir_val
      xdir = xdir_val
c .. get the field object id
      ISTAT = OMFfind_subobj(obj_id,
     .	OMFstr_to_name('out_fld'),
     .	OM_OBJ_RD, tmp_id)
      ISTAT = OMFget_obj_val(tmp_id,fld_id)
      RETURN
      END
      SUBROUTINE WRITE_MESSAGE(OBJ_ID,MESSAGE)
C***********************************************************************
C Write out a termination message, for the benefit of PC
C users who don't see stdout - connected to the UI
C***********************************************************************
      INCLUDE 'avs/omf.inc'
      INTEGER OBJ_ID(OIDSIZ),MESS_ID(OIDSIZ)
      CHARACTER*40 MESSAGE
      ISTAT = OMFfind_subobj(obj_id, OMFstr_to_name('message'),
     .				OM_OBJ_RD, MESS_ID)
      ISTAT = OMFset_str_val(MESS_ID,MESSAGE)
      RETURN
      END
      SUBROUTINE WRITE_TITLE(OBJ_ID,TITLE)
C***********************************************************************
C Write the title
C***********************************************************************
      INCLUDE 'avs/omf.inc'
      INTEGER OBJ_ID(OIDSIZ),TITLE_ID(OIDSIZ)
      CHARACTER*40 TITLE
      ISTAT = OMFfind_subobj(obj_id, OMFstr_to_name('title'),
     .				OM_OBJ_RD, title_id)
      ISTAT = OMFset_str_val(title_id,title)
      RETURN
      END
      SUBROUTINE WRITE_COORDS(FLD_ID,NMAP,COORD,NUMNOD)
C***********************************************************************
C Write the coordinates into the field structure
C***********************************************************************
      INCLUDE 'avs/omf.inc'
      INTEGER FLD_ID(OIDSIZ),NUMNOD,TMP_ID(OIDSIZ),NMAP(*)
      REAL COORD(3,*)
      REAL RBASE(1)
      INTEGER IBASE(1)
      ISTAT = OMFset_name_int_val(fld_id,
     .	OMFstr_to_name('nnodes'), NUMNOD)
      ISTAT = OMFset_name_int_val(fld_id,
     .	OMFstr_to_name('nspace'), 3)
      ISTAT = OMFfind_subobj(fld_id, OMFstr_to_name('coordinates'),
     .	OM_OBJ_RW, tmp_id)
      IPTR = OMFret_name_array_ptr(tmp_id,
     .	OMFstr_to_name('values'),
     .    OM_GET_ARRAY_WR, isize, itype)
      ioffs = ARRFoffset(IPTR,rbase,DTYPE_FLOAT)
c .. write the coordinates into the array
      DO 30 IN = 1,NUMNOD
        NOD = IN
        IF(IABS(NMAP(IN)) .LT. IN) NOD = -IN
        i = ioffs+(nod-1)*3
        rbase(i+1) = coord(1,NOD)
        rbase(i+2) = coord(2,NOD)
        rbase(i+3) = coord(3,NOD)
30	CONTINUE
      call ARRFfree(IPTR)
      RETURN
      END
      SUBROUTINE WRITE_ELEMENTS(FLD_ID,ELE,NUMELE)
C**********************************************************************
C Write the elements into the field structure
C**********************************************************************
      INCLUDE 'avs/omf.inc'
      INTEGER FLD_ID(OIDSIZ), TMP_ID(OIDSIZ), SET_ID(OIDSIZ),
     .	ELE(8,*), NUMELE
      INTEGER IBASE(1)
      ISTAT = OMFfind_subobj(FLD_ID,
     .    OMFstr_to_name('cell_set'),
     .    OM_OBJ_RD, set_id)
      ISTAT = OMFset_name_int_val(FLD_ID,
     .	OMFstr_to_name('ncell_sets'),1)
c .. init cell set 0
      ISTAT = OMFget_array_val(set_id,0,tmp_id,OM_OBJ_RW)
      ISTAT = OMFset_name_int_val(tmp_id,
     .	OMFstr_to_name('ncells'), NUMELE)
      ISTAT = OMFset_name_int_val(tmp_id,
     .	OMFstr_to_name('cell_ndim'), 3)
      ISTAT = OMFset_name_int_val(tmp_id,
     .	OMFstr_to_name('cell_order'), 1)
      ISTAT = OMFset_name_int_val(tmp_id,
     .	OMFstr_to_name('cell_nnodes'), 8)
      ISTAT = OMFset_name_int_val(tmp_id,
     .	OMFstr_to_name('poly_flag'), 0)
c .. get element pointer
      IPTR = OMFret_name_array_ptr(tmp_id,
     .	OMFstr_to_name('node_connect_list'),
     .	OM_GET_ARRAY_WR, ISIZE, ITYPE)
      ioffs = ARRFoffset(IPTR,ibase,DTYPE_INT)
c .. write the elements (zero-based) into the array
      DO 31 I = 1,NUMELE
      DO 31 J = 1,8
        KK = IOFFS+(I-1)*8+J
        IBASE(KK) = ELE(J,I)-1
31	CONTINUE
      call ARRFfree(IPTR)
      RETURN
      END
      SUBROUTINE INIT_FLDDATA(FLD_ID,NUMRES,CDA_ID,NDA_ID)
C**********************************************************************
C Initialise the cell data structure - set numres and return IDs for
C cell and node data
C**********************************************************************
      INCLUDE 'avs/omf.inc'
      INTEGER FLD_ID(OIDSIZ), NUMRES, CDA_ID(OIDSIZ),
     .    NDA_ID(OIDSIZ), TMP_ID(OIDSIZ), SET_ID(OIDSIZ),
     .    VAL_ID(OIDSIZ)
      INTEGER ISTAT
c .. get the cell set
      ISTAT = OMFfind_subobj(fld_id,
     .	OMFstr_to_name('cell_set'), OM_OBJ_RD, set_id)
c .. set ncell_data
      ISTAT = OMFget_array_val(set_id,0,val_id,OM_OBJ_RW)
      ISTAT = OMFset_name_int_val(val_id,
     .	OMFstr_to_name('ncell_data'),
     .	NUMRES+1)
c .. get the cell data array id
      ISTAT = OMFfind_subobj(val_id,
     .	OMFstr_to_name('cell_data'),
     .	OM_OBJ_RD, cda_id)
c .. set nnode_data
      ISTAT = OMFset_name_int_val(fld_id,
     .	OMFstr_to_name('nnode_data'), NUMRES)
c .. get the node data array id
      ISTAT = OMFfind_subobj(fld_id,
     .	OMFstr_to_name('node_data'),
     .	OM_OBJ_WR, nda_id)
c ..
      RETURN
      END
      SUBROUTINE WRITE_DATA(TMP_ID,INDEX,NAME,NVALS,VALS)
C**********************************************************************
C Write an array into a data object
C**********************************************************************
      INCLUDE 'avs/omf.inc'
      INTEGER TMP_ID(OIDSIZ),INDEX,VAL_ID(OIDSIZ),IPTR,NVALS
	INTEGER	NAME_ID(OIDSIZ)
      REAL RBASE(1),VALS(NVALS)
      CHARACTER*(*) NAME
c .. set array
      ISTAT = OMFget_array_val(TMP_ID,INDEX,VAL_ID,OM_OBJ_RW)
      IPTR = OMFret_name_array_ptr(VAL_ID,
     .	OMFstr_to_name('values'),
     .	OM_GET_ARRAY_WR,
     .	ISIZE, ITYPE)
      IOFFS = ARRFoffset(IPTR,RBASE,DTYPE_FLOAT)
      DO 33 I = 1,NVALS
      	RBASE(IOFFS+I) = VALS(I)
33	CONTINUE
      CALL ARRFfree(IPTR)
c .. set data label
      ISTAT = OMFset_name_str_val(VAL_ID,
     .    OMFstr_to_name('labels'), NAME)
      RETURN
      END
      SUBROUTINE FCLSNP(IOPHI,IOPHDA,IONEP,IOPAT,IOGRD)
C***********************************************************************
C  Author: CYC
C  Creation Date: 18th July 1994
C  Amendment History
c .. Nature      Date      Author
C
C  Function of Subroutine:  Closes the files.
C
C  Called by: MAINNP
C
C  Limitations
C
C  Source of further information
C
C  Examples of use
C
C  Input Parameters
C
C  Output Parameters
C
C***********************************************************************
C
      CLOSE(IOPHI)
      CLOSE(IOPHDA)
      CLOSE(IONEP)
      IF (IOPAT.NE.0) THEN
      CLOSE(IOPAT)
	ENDIF
      CLOSE(IOGRD)

      RETURN
      END
      SUBROUTINE WELENP(IONEP,IELEM,NUMELE)
C***********************************************************************
C  Author: CYC
C  Creation Date: 18th July 1994
C  Amendment History
c .. Nature      Date      Author
C
C  Function of Subroutine:  Write out the element definition.
C
C  Called by: MAINNP
C
C  Limitations
C
C  Source of further information
C
C  Examples of use
C
C  Input Parameters
C
C  Output Parameters
C
C***********************************************************************
C
C
      INTEGER IELEM(8,*)
C
      WRITE(IONEP,*)NUMELE
      DO 100 IE = 1,NUMELE
        WRITE(IONEP,'(9(1X,I5))')IE,(IELEM(IN,IE),IN=1,8)
 100  CONTINUE
C
      RETURN
      END
      SUBROUTINE WNODNP(IONEP,NMAP,COORD,NUMNOD)
C***********************************************************************
C  Author: CYC
C  Creation Date: 18th July 1994
C  Amendment History
c .. Nature      Date      Author
C
C  Function of Subroutine:  Converts the PHOENICS PHI(DA),PATGEO and
c ..                       XYZ(DA) files into the NEP file.
C
C  Called by: NEPTRN
C
C  Limitations
C
C  Source of further information
C
C  Examples of use
C
C  Input Parameters
C
C  Output Parameters
C
C***********************************************************************
      REAL COORD(3,*)
      INTEGER NMAP(*)
C
      WRITE(IONEP,*)NUMNOD
      DO 100 IN = 1,NUMNOD
        NOD = IN
        IF(IABS(NMAP(IN)) .LT. IN) NOD = -IN
c ..    IF(NMAP(IN) .GT. 0 .AND. NMAP(IN) .NE. IN) NOD = -IN
        WRITE(IONEP,*)NOD,COORD(1,IN),COORD(2,IN),COORD(3,IN)
 100  CONTINUE
C
      RETURN
      END
      SUBROUTINE OPHINP(IONEP,IOPHI,IOPHDA,DIRECT,LBYTE,INT,PHI,LENREF,
     +                    BFC,CARTES,NUMRES,LOG,NAME,IERR)
C***********************************************************************
C
C  Author: CYC
C  Creation Date: 18th July 1994
C  Amendment History
c .. Nature      Date      Author
C
C  Function of Subroutine:  Open the phi file. If ASCII convert to DA
C
C  Called by: MAINPN
C
C  Limitations
C
C  Source of further information
C
C  Examples of use
C
C  Input Parameters
c .. IONEP  - NEP file unit number
c .. IOPHI  - PHI file unit number
c .. IOPHDA - PHIDA unit number
c .. DIRECT - PHI file type
c .. LBYTE  - Number of bytes / words
c .. INT     - Work array for reading Integers
c .. PHI     - Work array for reading Reals
c .. LENREF  - Number of reals / record
C
C  Output Parameters
c .. BFC    - Grid is BFC flag
c .. CARTES  - Polar or Cartesian flag
c .. NUMRES  - Number of results found
c .. LOG     - Array of pointers to results found
c .. NAME    - Name of results found
c .. IERR    - Error flag
C
C***********************************************************************
      EXTERNAL ERRNP
C
      REAL REC(6)
      REAL PHI(LENREF)
      LOGICAL LOG(50),LOG1,LOG2,LOG3,LOG4
      CHARACTER*40 TITLE
      CHARACTER*4 NAME(50)
      INTEGER INT(LENREF)
C
      LOGICAL DIRECT,BFC,CARTES
C
      LEN = LBYTE*LENREF
      IF(DIRECT) THEN
c ..   Open with correct record length
        OPEN(IOPHDA,FILE='PHIDA',STATUS='OLD',RECL=LEN,
     +        ACCESS='DIRECT',ERR=900)
C
c ..   ASCII file
      ELSE
c ..   Open the PHI file
        OPEN(IOPHI,FILE='PHI',STATUS='OLD',ERR=800)
c ..   Convert this to a direct access file
        OPEN(IOPHDA,STATUS='SCRATCH',RECL=LEN,
     +       ACCESS='DIRECT',ERR=900)
C
        READ (IOPHI,1,ERR=700)TITLE
 1      FORMAT(A40)
        WRITE(IOPHDA,REC=1)TITLE
C
        READ(IOPHI,2,ERR=700)LOG1,LOG2,LOG3,LOG4
 2      FORMAT(1X,4L1)
        WRITE(IOPHDA,REC=2)LOG1,LOG2,LOG3,LOG4
C
        READ(IOPHI,3,ERR=700)(INT(J),J=1,11)
 3      FORMAT(1X,8I10,/,1X,4I10)
        WRITE(IOPHDA,REC=3)INT
        NX = INT(1)
        NY = INT(2)
        NZ = INT(3)
C
        READ(IOPHI,4,ERR=700)PHI(1)
 4      FORMAT(E13.2)
        WRITE(IOPHDA,REC=4)PHI(1)
C
        READ(IOPHI,5,ERR=700)NAME
 5      FORMAT(1X,19A4,/,1X,19A4,/,1X,12A4)
        WRITE(IOPHDA,REC=5)NAME
C
        READ(IOPHI,6,ERR=700)(PHI(I),I=1,NX)
 6      FORMAT(43(6(E13.0:),/))
        WRITE(IOPHDA,REC=6)(PHI(I),I=1,NX)
C
        READ(IOPHI,7,ERR=700)(PHI(I),I=1,NY)
 7      FORMAT(43(6(E13.0:),/))
        WRITE(IOPHDA,REC=7)(PHI(I),I=1,NY)
C
        READ(IOPHI,8,ERR=700)(PHI(I),I=1,NZ)
 8      FORMAT(43(6(E13.0:),/))
        WRITE(IOPHDA,REC=8)(PHI(I),I=1,NZ)
C
        READ(IOPHI,9,ERR=700)(PHI(I),I=1,NZ)
 9      FORMAT(43(6(E13.0:),/))
        WRITE(IOPHDA,REC=9)(PHI(I),I=1,NZ)
C
        READ(IOPHI,10,ERR=700)LOG
 10     FORMAT(1X,50L1)
        WRITE(IOPHDA,REC=10)LOG
C
C This is a do while loop
C
        IREC = 10
        NLAYER = NX*NY
        NUMREC = NLAYER/6
        NREST  = MOD(NLAYER,6)
 100    CONTINUE
        IPHI = 1
c ..   For All the ascii records but one for this layer
        DO 110 NASC = 1,NUMREC
          READ(IOPHI,115,END=200,ERR=700)REC
 115      FORMAT(6E13.0)
          DO 120 IN = 1,6
            PHI(IPHI) = REC(IN)
c ..       Full binary record write it out
            IF (IPHI .LT. 256) THEN
              IPHI = IPHI + 1
            ELSE
              IREC = IREC + 1
              WRITE(IOPHDA,REC=IREC)PHI
              IPHI = 1
            ENDIF
 120      CONTINUE
 110    CONTINUE
c ..
c ..   Last record for the layer
        IF (NREST .NE. 0) THEN
          READ(IOPHI,125,END=200,ERR=700)(REC(I),I=1,NREST)
 125      FORMAT(6E13.0:)
          DO 130 IN = 1,NREST
            PHI(IPHI) = REC(IN)
            IF (IPHI .LT. 256) THEN
              IPHI = IPHI + 1
            ELSE
              IREC = IREC + 1
              WRITE(IOPHDA,REC=IREC)PHI
              IPHI = 1
            ENDIF
 130      CONTINUE
        ENDIF
c ..   Write out incomplete record
        IF (IPHI .NE. 1)THEN
          IREC = IREC + 1
          WRITE(IOPHDA,REC=IREC)(PHI(I),I=1,(IPHI-1))
        ENDIF
C
        GOTO 100
 200    CONTINUE
        CLOSE(IOPHI)
      ENDIF
      READ(IOPHDA,REC=1,ERR=750)TITLE
      READ(IOPHDA,REC=2,ERR=750)CARTES,LOG2,BFC,LOG4
      READ(IOPHDA,REC=5,ERR=750)NAME
      READ(IOPHDA,REC=10,ERR=750)LOG
      GOTO 999
 700  CONTINUE
      IERR = 5
      CALL ERRNP(IONEP,IERR,' ')
      GOTO 999
 750  CONTINUE
      IERR = 6
      CALL ERRNP(IONEP,IERR,' ')
      GOTO 999
 800  CONTINUE
      IERR = 3
      CALL ERRNP(IONEP,IERR,' ')
      GOTO 999
 900  CONTINUE
      IERR = 4
      CALL ERRNP(IONEP,IERR,' ')
 999  CONTINUE
      RETURN
      END

      SUBROUTINE OPATNP(IONEP,IOPAT,PATNAM,IERR)
C***********************************************************************
C  Author: CYC
C  Creation Date: 18th July 1994
C  Amendment History
c .. Nature      Date      Author
C
C  Function of Subroutine:  Open the PATGEO file
C
C  Called by: NEPTRN
C
C  Limitations
C
C  Source of further information
C
C  Examples of use
C
C  Input Parameters
c ..  IONEP - NEP File Unit number
c ..  IOPAT - PATGEO File Unit number
C
C  Output Parameters
c ..  IERR  - File opening error flag
C
C***********************************************************************
      CHARACTER*(*) PATNAM
      IERR= 0
      OPEN(UNIT=IOPAT,FILE=PATNAM,STATUS='OLD',ERR=10)
      GOTO 20
 10   IERR = 2
	IOPAT = 0
      CALL ERRNP(IONEP,IERR,' ')
 20   CONTINUE
      RETURN
      END
      SUBROUTINE WRTLNK(IONEP,LINK,NLINK,IJKBLK,NBLOCK)
C***********************************************************************
C  Author: CYC
C  Creation Date: 18th July 1994
C  Amendment History
c .. Nature      Date      Author
C
C  Function of Subroutine:  Writes links in local IJK to end of NEP file.
C
C  Called by: NEPTRN
C
C  Limitations
C
C  Source of further information
C
C  Examples of use
C
C  Input Parameters
C    IONEP  - NEP file unit number
C    LINK   - Link definitions
C    NLINK  - Number of links
C    IJKBLK - Block limits
C    NBLOCK - Number of blocks
C
C  Output Parameters
C
C***********************************************************************
      INTEGER IJKBLK(6,NBLOCK)
      INTEGER LINK(11,NLINK)
C
c .. For each link
c ..
      WRITE(IONEP,*)NLINK
      DO 100 IL = 1,NLINK
c ..   Find element offset for first block in link
        IB1 = IABS(LINK(1,IL))
        IEOF = IEOFF(IB1,IJKBLK)
c ..   Write link definition
        IB2 = IABS(LINK(2,IL))
        ID1 = LINK(3,IL)
        ID1S = LINK(4,IL) - IJKBLK(ID1*2-1,IB1) + 1
        ID1F = LINK(5,IL) - IJKBLK(ID1*2-1,IB1) + 1
        ID2 = LINK(6,IL)
        ID2S = LINK(7,IL) - IJKBLK(ID2*2-1,IB1) + 1
        ID2F = LINK(8,IL) - IJKBLK(ID2*2-1,IB1) + 1
        ID3 = LINK(9,IL)
        ID3S = LINK(10,IL) - IJKBLK(ID3*2-1,IB1) + 1
        ID3F = LINK(11,IL) - IJKBLK(ID3*2-1,IB1) + 1
        WRITE(IONEP,10)IB1,IB2,ID1,ID1S,ID1F,ID2,ID2S,ID2F,ID3,ID3S,ID3F
 10     FORMAT(2(I8,1X),3(I1,1X,I8,1X,I8,1X))
C
 100  CONTINUE
      RETURN
      END

ext: f

max_stars_repo_name: avs/express-community

max_issues_repo_name: avs/express-community

max_forks_repo_name: avs/express-community


text:
subroutine MultMassMatrix(a,array,idofn,ndofn)
   use typre
   use Mod_Memor
   use Mod_ParallelSystemInterface
   use Mod_Element
   use Mod_Mesh
   implicit none
   class(ParallelSystemInterface),pointer :: LinearSystem => NULL()
   class(FemMesh) :: a
   class(FiniteElement), pointer :: e => NULL()
   integer(ip) :: inode,jnode,igaus,ielem,idofn,ipoin,npoin,nelem,npoinLocal,ndofn
   real(rp)    :: array(ndofn,a%npoin),adiag,dvol,vecout(a%npoin),vecin(a%npoin)
   real(rp), allocatable :: elmat(:,:),elrhs(:)
   character(150) :: auxstring, auxstring2

   call a%ElementAlloc(e,a%Memor,'DefaultRule','L2Projector')
   call a%Memor%Alloc(e%mnode,e%mnode,elmat,'elmat','L2Projector')
   call a%Memor%Alloc(e%mnode,elrhs,'elrhs','L2Projector')
   call a%GetNelem(nelem)
   call a%GetNpoinLocal(npoinLocal)

   auxstring = trim(a%InputFolder)//'/'//adjustl(trim(a%namda))//'.dom.sol '
   auxstring2 = 'dom '

   call a%ParallelLibrary%CreateSystem(LinearSystem, a%Memor)
   call a%InitializeSystem(1_ip,0_ip,LinearSystem,a%Memor,auxstring,auxstring2,a%lun_solve_dom)

   call LinearSystem%ToZero
   vecin(:) = array(idofn,:)

   !Element loop
   elements : do ielem = 1,nelem
      !Load Element
      call a%ElementLoad(ielem,e)
      call e%elmdel

      elmat = 0.0_rp

      gauss_points : do igaus=1,e%pgaus
         e%igaus = igaus
         call e%elmder
         dvol = e%weigp(e%igaus)*e%detjm
         do inode = 1,e%pnode
            do jnode = 1,e%pnode
               elmat(inode,jnode) = elmat(inode,jnode) + dvol*e%shape(inode,e%igaus)*e%shape(jnode,e%igaus)
            enddo
         enddo
       enddo gauss_points

       !Assembly
       call LinearSystem%Assembly(e,elmat,elrhs)
   enddo elements

   call LinearSystem%Matvec(vecin,vecout)
   array(idofn,:) = vecout

   !Communicate Ghosts
   call a%ArrayCommunicator%GhostCommunicate(1,array(idofn,:))

   call a%Memor%Dealloc(e%mnode,e%mnode,elmat,'elmat','L2Projector')
   call a%Memor%Dealloc(e%mnode,elrhs,'elrhs','L2Projector')
   call a%ElementDealloc(e,a%Memor,'DefaultRule','Projector')
   call LinearSystem%Deallocate
   call a%ParallelLibrary%DeallocateSystem(LinearSystem,a%Memor)
end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
C Copyright 2021 Dennis Decker Jensen
C Date: 9 May 2021
C Purpose: Prime number test
C Tectonics: gfortran --std=f95 -ffixed-form -c primetest.f
      LOGICAL FUNCTION PRIME(NUM)
         IROOT = INT(SQRT(REAL(NUM)))

         DO 10 I=2,IROOT
            IF (MOD(NUM,I).EQ.0) THEN
               PRIME=.FALSE.
               RETURN
            ENDIF
   10    CONTINUE
         PRIME=.TRUE.
      END

max_stars_repo_name: dennisdjensen/Sketchbook

max_issues_repo_name: dennisdjensen/sketchbook

max_forks_repo_name: dennisdjensen/sketchbook


text:
Talbott Solar in Davis, California, provides residential, small business, and rental property solar power solar electric and heating systems.

Talbott Solar was responsible for the solar panels at the Whole Earth Festivals Main Stage.  An exposition at the Explorit Science Center showcasing a doll house with solar panels was funded by Talbott Solar.

Talbott Solar also installs radiant heating and cooling celling systems, solar hot water heaters, tankless hot water heaters, heat pump water heaters, and hydronic hot water heaters.

Talbott Solar is owned and operated by Dean Newberry. A great http://search.davisenterprise.com/display.php?id66531 article about his life and his newly invented with a patent pending radiant pannels was recently in the Davis Enterprise.

20091123 17:05:33 nbsp Installed SunPower systems, very good looking, but expensive. Users/johndeloro

20100826 08:27:32 nbsp Great products! Originally installed tankless water heater was undersized for my house. Received NO assistance from owner and sales rep in trouble shooting the defective pool pump controller that was installed by Talbott  disappointing after sales service. Users/ReneTsang

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 hecmw_solver_scaling_66
  use hecmw_util
  use m_hecmw_comm_f
  use hecmw_matrix_misc
  implicit none

  private
  real(kind=kreal), private, allocatable :: scale(:)

  public :: hecmw_solver_scaling_fw_66
  public :: hecmw_solver_scaling_bk_66

contains

  subroutine hecmw_solver_scaling_fw_66(hecMESH, hecMAT, COMMtime)
    implicit none
    type (hecmwST_local_mesh), intent(in) :: hecMESH
    type (hecmwST_matrix), intent(inout) :: hecMAT
    real(kind=kreal), intent(inout) :: COMMtime
    integer(kind=kint) :: N, NP, NDOF
    real(kind=kreal), pointer :: D(:), AL(:), AU(:), B(:)
    integer(kind=kint), pointer :: INL(:), IAL(:), INU(:), IAU(:)
    integer(kind=kint) :: i, k, ip1, ip2, ip3, ip4, ip5, ip6
    integer(kind=kint) :: iq1, iq2, iq3, iq4, iq5, iq6
    integer(kind=kint) :: isL, ieL, isU, ieU, inod
    real(kind=kreal) :: START_TIME, END_TIME

    if (hecmw_mat_get_scaling(hecMAT).eq.0) return

    N = hecMAT%N
    NP = hecMAT%NP
    NDOF = hecMAT%NDOF
    D => hecMAT%D
    AL => hecMAT%AL
    AU => hecMAT%AU
    INL => hecMAT%indexL
    IAL => hecMAT%itemL
    INU => hecMAT%indexU
    IAU => hecMAT%itemU
    B => hecMAT%B

    allocate(scale(NDOF*NP))

    do i= 1, N
      scale (6*i-5)= 1.d0/dsqrt(dabs(D(36*i-35)))
      scale (6*i-4)= 1.d0/dsqrt(dabs(D(36*i-28)))
      scale (6*i-3)= 1.d0/dsqrt(dabs(D(36*i-21)))
      scale (6*i-2)= 1.d0/dsqrt(dabs(D(36*i-14)))
      scale (6*i-1)= 1.d0/dsqrt(dabs(D(36*i-7)))
      scale (6*i  )= 1.d0/dsqrt(dabs(D(36*i  )))
    enddo

    START_TIME= HECMW_WTIME()
    call hecmw_update_3_R (hecMESH, scale, hecMESH%n_node)
    END_TIME= HECMW_WTIME()
    COMMtime = COMMtime + END_TIME - START_TIME

    do i= 1, NP
      ip1= 6*i-5
      ip2= 6*i-4
      ip3= 6*i-3
      ip4= 6*i-2
      ip5= 6*i-1
      ip6= 6*i

      D(36*i-35)= D(36*i-35)*scale(ip1)*scale(ip1)
      D(36*i-34)= D(36*i-34)*scale(ip1)*scale(ip2)
      D(36*i-33)= D(36*i-33)*scale(ip1)*scale(ip3)
      D(36*i-32)= D(36*i-32)*scale(ip1)*scale(ip4)
      D(36*i-31)= D(36*i-31)*scale(ip1)*scale(ip5)
      D(36*i-30)= D(36*i-30)*scale(ip1)*scale(ip6)

      D(36*i-29)= D(36*i-29)*scale(ip2)*scale(ip1)
      D(36*i-28)= D(36*i-28)*scale(ip2)*scale(ip2)
      D(36*i-27)= D(36*i-27)*scale(ip2)*scale(ip3)
      D(36*i-26)= D(36*i-26)*scale(ip2)*scale(ip4)
      D(36*i-25)= D(36*i-25)*scale(ip2)*scale(ip5)
      D(36*i-24)= D(36*i-24)*scale(ip2)*scale(ip6)

      D(36*i-23)= D(36*i-23)*scale(ip3)*scale(ip1)
      D(36*i-22)= D(36*i-22)*scale(ip3)*scale(ip2)
      D(36*i-21)= D(36*i-21)*scale(ip3)*scale(ip3)
      D(36*i-20)= D(36*i-20)*scale(ip3)*scale(ip4)
      D(36*i-19)= D(36*i-19)*scale(ip3)*scale(ip5)
      D(36*i-18)= D(36*i-18)*scale(ip3)*scale(ip6)

      D(36*i-17)= D(36*i-17)*scale(ip4)*scale(ip1)
      D(36*i-16)= D(36*i-16)*scale(ip4)*scale(ip2)
      D(36*i-15)= D(36*i-15)*scale(ip4)*scale(ip3)
      D(36*i-14)= D(36*i-14)*scale(ip4)*scale(ip4)
      D(36*i-13)= D(36*i-13)*scale(ip4)*scale(ip5)
      D(36*i-12)= D(36*i-12)*scale(ip4)*scale(ip6)

      D(36*i-11)= D(36*i-11)*scale(ip5)*scale(ip1)
      D(36*i-10)= D(36*i-10)*scale(ip5)*scale(ip2)
      D(36*i-9 )= D(36*i-9 )*scale(ip5)*scale(ip3)
      D(36*i-8 )= D(36*i-8 )*scale(ip5)*scale(ip4)
      D(36*i-7 )= D(36*i-7 )*scale(ip5)*scale(ip5)
      D(36*i-6 )= D(36*i-6 )*scale(ip5)*scale(ip6)

      D(36*i-5 )= D(36*i-5 )*scale(ip6)*scale(ip1)
      D(36*i-4 )= D(36*i-4 )*scale(ip6)*scale(ip2)
      D(36*i-3 )= D(36*i-3 )*scale(ip6)*scale(ip3)
      D(36*i-2 )= D(36*i-2 )*scale(ip6)*scale(ip4)
      D(36*i-1 )= D(36*i-1 )*scale(ip6)*scale(ip5)
      D(36*i   )= D(36*i   )*scale(ip6)*scale(ip6)

      isL= INL(i-1) + 1
      ieL= INL(i  )
      !*voption indep (IAL,AL,SCALE)
      do k= isL, ieL
        inod= IAL(k)
        iq1= 6*inod - 5
        iq2= 6*inod - 4
        iq3= 6*inod - 3
        iq4= 6*inod - 2
        iq5= 6*inod - 1
        iq6= 6*inod

        AL(36*k-35)= AL(36*k-35)*scale(ip1)*scale(iq1)
        AL(36*k-34)= AL(36*k-34)*scale(ip1)*scale(iq2)
        AL(36*k-33)= AL(36*k-33)*scale(ip1)*scale(iq3)
        AL(36*k-32)= AL(36*k-32)*scale(ip1)*scale(iq4)
        AL(36*k-31)= AL(36*k-31)*scale(ip1)*scale(iq5)
        AL(36*k-30)= AL(36*k-30)*scale(ip1)*scale(iq6)

        AL(36*k-29)= AL(36*k-29)*scale(ip2)*scale(iq1)
        AL(36*k-28)= AL(36*k-28)*scale(ip2)*scale(iq2)
        AL(36*k-27)= AL(36*k-27)*scale(ip2)*scale(iq3)
        AL(36*k-26)= AL(36*k-26)*scale(ip2)*scale(iq4)
        AL(36*k-25)= AL(36*k-25)*scale(ip2)*scale(iq5)
        AL(36*k-24)= AL(36*k-24)*scale(ip2)*scale(iq6)

        AL(36*k-23)= AL(36*k-23)*scale(ip3)*scale(iq1)
        AL(36*k-22)= AL(36*k-22)*scale(ip3)*scale(iq2)
        AL(36*k-21)= AL(36*k-21)*scale(ip3)*scale(iq3)
        AL(36*k-20)= AL(36*k-20)*scale(ip3)*scale(iq4)
        AL(36*k-19)= AL(36*k-19)*scale(ip3)*scale(iq5)
        AL(36*k-18)= AL(36*k-18)*scale(ip3)*scale(iq6)

        AL(36*k-17)= AL(36*k-17)*scale(ip4)*scale(iq1)
        AL(36*k-16)= AL(36*k-16)*scale(ip4)*scale(iq2)
        AL(36*k-15)= AL(36*k-15)*scale(ip4)*scale(iq3)
        AL(36*k-14)= AL(36*k-14)*scale(ip4)*scale(iq4)
        AL(36*k-13)= AL(36*k-13)*scale(ip4)*scale(iq5)
        AL(36*k-12)= AL(36*k-12)*scale(ip4)*scale(iq6)

        AL(36*k-11)= AL(36*k-11)*scale(ip5)*scale(iq1)
        AL(36*k-10)= AL(36*k-10)*scale(ip5)*scale(iq2)
        AL(36*k-9 )= AL(36*k-9 )*scale(ip5)*scale(iq3)
        AL(36*k-8 )= AL(36*k-8 )*scale(ip5)*scale(iq4)
        AL(36*k-7 )= AL(36*k-7 )*scale(ip5)*scale(iq5)
        AL(36*k-6 )= AL(36*k-6 )*scale(ip5)*scale(iq6)

        AL(36*k-5 )= AL(36*k-5 )*scale(ip6)*scale(iq1)
        AL(36*k-4 )= AL(36*k-4 )*scale(ip6)*scale(iq2)
        AL(36*k-3 )= AL(36*k-3 )*scale(ip6)*scale(iq3)
        AL(36*k-2 )= AL(36*k-2 )*scale(ip6)*scale(iq4)
        AL(36*k-1 )= AL(36*k-1 )*scale(ip6)*scale(iq5)
        AL(36*k   )= AL(36*k   )*scale(ip6)*scale(iq6)
      enddo

      isU= INU(i-1) + 1
      ieU= INU(i  )
      !*voption indep (IAU,AU,SCALE)
      do k= isU, ieU
        inod= IAU(k)
        iq1= 6*inod - 5
        iq2= 6*inod - 4
        iq3= 6*inod - 3
        iq4= 6*inod - 2
        iq5= 6*inod - 1
        iq6= 6*inod

        AU(36*k-35)= AU(36*k-35)*scale(ip1)*scale(iq1)
        AU(36*k-34)= AU(36*k-34)*scale(ip1)*scale(iq2)
        AU(36*k-33)= AU(36*k-33)*scale(ip1)*scale(iq3)
        AU(36*k-32)= AU(36*k-32)*scale(ip1)*scale(iq4)
        AU(36*k-31)= AU(36*k-31)*scale(ip1)*scale(iq5)
        AU(36*k-30)= AU(36*k-30)*scale(ip1)*scale(iq6)

        AU(36*k-29)= AU(36*k-29)*scale(ip2)*scale(iq1)
        AU(36*k-28)= AU(36*k-28)*scale(ip2)*scale(iq2)
        AU(36*k-27)= AU(36*k-27)*scale(ip2)*scale(iq3)
        AU(36*k-26)= AU(36*k-26)*scale(ip2)*scale(iq4)
        AU(36*k-25)= AU(36*k-25)*scale(ip2)*scale(iq5)
        AU(36*k-24)= AU(36*k-24)*scale(ip2)*scale(iq6)

        AU(36*k-23)= AU(36*k-23)*scale(ip3)*scale(iq1)
        AU(36*k-22)= AU(36*k-22)*scale(ip3)*scale(iq2)
        AU(36*k-21)= AU(36*k-21)*scale(ip3)*scale(iq3)
        AU(36*k-20)= AU(36*k-20)*scale(ip3)*scale(iq4)
        AU(36*k-19)= AU(36*k-19)*scale(ip3)*scale(iq5)
        AU(36*k-18)= AU(36*k-18)*scale(ip3)*scale(iq6)

        AU(36*k-17)= AU(36*k-17)*scale(ip4)*scale(iq1)
        AU(36*k-16)= AU(36*k-16)*scale(ip4)*scale(iq2)
        AU(36*k-15)= AU(36*k-15)*scale(ip4)*scale(iq3)
        AU(36*k-14)= AU(36*k-14)*scale(ip4)*scale(iq4)
        AU(36*k-13)= AU(36*k-13)*scale(ip4)*scale(iq5)
        AU(36*k-12)= AU(36*k-12)*scale(ip4)*scale(iq6)

        AU(36*k-11)= AU(36*k-11)*scale(ip5)*scale(iq1)
        AU(36*k-10)= AU(36*k-10)*scale(ip5)*scale(iq2)
        AU(36*k-9 )= AU(36*k-9 )*scale(ip5)*scale(iq3)
        AU(36*k-8 )= AU(36*k-8 )*scale(ip5)*scale(iq4)
        AU(36*k-7 )= AU(36*k-7 )*scale(ip5)*scale(iq5)
        AU(36*k-6 )= AU(36*k-6 )*scale(ip5)*scale(iq6)

        AU(36*k-5 )= AU(36*k-5 )*scale(ip6)*scale(iq1)
        AU(36*k-4 )= AU(36*k-4 )*scale(ip6)*scale(iq2)
        AU(36*k-3 )= AU(36*k-3 )*scale(ip6)*scale(iq3)
        AU(36*k-2 )= AU(36*k-2 )*scale(ip6)*scale(iq4)
        AU(36*k-1 )= AU(36*k-1 )*scale(ip6)*scale(iq5)
        AU(36*k   )= AU(36*k   )*scale(ip6)*scale(iq6)
      enddo
    enddo
    !*voption indep (B,SCALE)
    do i= 1, N
      B(6*i-5)= B(6*i-5) * scale(6*i-5)
      B(6*i-4)= B(6*i-4) * scale(6*i-4)
      B(6*i-3)= B(6*i-3) * scale(6*i-3)
      B(6*i-2)= B(6*i-2) * scale(6*i-2)
      B(6*i-1)= B(6*i-1) * scale(6*i-1)
      B(6*i  )= B(6*i  ) * scale(6*i  )
    enddo
  end subroutine hecmw_solver_scaling_fw_66

  subroutine hecmw_solver_scaling_bk_66(hecMAT)
    use hecmw_util
    implicit none
    type (hecmwST_matrix), intent(inout) :: hecMAT
    integer(kind=kint) :: N, NP, NDOF
    real(kind=kreal), pointer :: D(:), AL(:), AU(:), B(:), X(:)
    integer(kind=kint), pointer :: INL(:), IAL(:), INU(:), IAU(:)
    integer(kind=kint) :: i, k, ip1, ip2, ip3, ip4, ip5, ip6
    integer(kind=kint) :: iq1, iq2, iq3, iq4, iq5, iq6
    integer(kind=kint) :: isL, ieL, isU, ieU, inod

    if (hecmw_mat_get_scaling(hecMAT).eq.0) return

    N = hecMAT%N
    NP = hecMAT%NP
    NDOF = hecMAT%NDOF
    D => hecMAT%D
    AL => hecMAT%AL
    AU => hecMAT%AU
    INL => hecMAT%indexL
    IAL => hecMAT%itemL
    INU => hecMAT%indexU
    IAU => hecMAT%itemU
    B => hecMAT%B
    X => hecMAT%X

    !*voption indep (X,B,SCALE)
    do i= 1, N
      X(6*i-5)= X(6*i-5) * scale(6*i-5)
      X(6*i-4)= X(6*i-4) * scale(6*i-4)
      X(6*i-3)= X(6*i-3) * scale(6*i-3)
      X(6*i-2)= X(6*i-2) * scale(6*i-2)
      X(6*i-1)= X(6*i-1) * scale(6*i-1)
      X(6*i  )= X(6*i  ) * scale(6*i  )

      B(6*i-5)= B(6*i-5) / scale(6*i-5)
      B(6*i-4)= B(6*i-4) / scale(6*i-4)
      B(6*i-3)= B(6*i-3) / scale(6*i-3)
      B(6*i-2)= B(6*i-2) / scale(6*i-2)
      B(6*i-1)= B(6*i-1) / scale(6*i-1)
      B(6*i  )= B(6*i  ) / scale(6*i  )
    enddo

    do i= 1, NP
      ip1= 6*i-5
      ip2= 6*i-4
      ip3= 6*i-3
      ip4= 6*i-2
      ip5= 6*i-1
      ip6= 6*i

      D(36*i-35)= D(36*i-35)/(scale(ip1)*scale(ip1))
      D(36*i-34)= D(36*i-34)/(scale(ip1)*scale(ip2))
      D(36*i-33)= D(36*i-33)/(scale(ip1)*scale(ip3))
      D(36*i-32)= D(36*i-32)/(scale(ip1)*scale(ip4))
      D(36*i-31)= D(36*i-31)/(scale(ip1)*scale(ip5))
      D(36*i-30)= D(36*i-30)/(scale(ip1)*scale(ip6))

      D(36*i-29)= D(36*i-29)/(scale(ip2)*scale(ip1))
      D(36*i-28)= D(36*i-28)/(scale(ip2)*scale(ip2))
      D(36*i-27)= D(36*i-27)/(scale(ip2)*scale(ip3))
      D(36*i-26)= D(36*i-26)/(scale(ip2)*scale(ip4))
      D(36*i-25)= D(36*i-25)/(scale(ip2)*scale(ip5))
      D(36*i-24)= D(36*i-24)/(scale(ip2)*scale(ip6))

      D(36*i-23)= D(36*i-23)/(scale(ip3)*scale(ip1))
      D(36*i-22)= D(36*i-22)/(scale(ip3)*scale(ip2))
      D(36*i-21)= D(36*i-21)/(scale(ip3)*scale(ip3))
      D(36*i-20)= D(36*i-20)/(scale(ip3)*scale(ip4))
      D(36*i-19)= D(36*i-19)/(scale(ip3)*scale(ip5))
      D(36*i-18)= D(36*i-18)/(scale(ip3)*scale(ip6))

      D(36*i-17)= D(36*i-17)/(scale(ip4)*scale(ip1))
      D(36*i-16)= D(36*i-16)/(scale(ip4)*scale(ip2))
      D(36*i-15)= D(36*i-15)/(scale(ip4)*scale(ip3))
      D(36*i-14)= D(36*i-14)/(scale(ip4)*scale(ip4))
      D(36*i-13)= D(36*i-13)/(scale(ip4)*scale(ip5))
      D(36*i-12)= D(36*i-12)/(scale(ip4)*scale(ip6))

      D(36*i-11)= D(36*i-11)/(scale(ip5)*scale(ip1))
      D(36*i-10)= D(36*i-10)/(scale(ip5)*scale(ip2))
      D(36*i-9 )= D(36*i-9 )/(scale(ip5)*scale(ip3))
      D(36*i-8 )= D(36*i-8 )/(scale(ip5)*scale(ip4))
      D(36*i-7 )= D(36*i-7 )/(scale(ip5)*scale(ip5))
      D(36*i-6 )= D(36*i-6 )/(scale(ip5)*scale(ip6))

      D(36*i-5 )= D(36*i-5 )/(scale(ip6)*scale(ip1))
      D(36*i-4 )= D(36*i-4 )/(scale(ip6)*scale(ip2))
      D(36*i-3 )= D(36*i-3 )/(scale(ip6)*scale(ip3))
      D(36*i-2 )= D(36*i-2 )/(scale(ip6)*scale(ip4))
      D(36*i-1 )= D(36*i-1 )/(scale(ip6)*scale(ip5))
      D(36*i   )= D(36*i   )/(scale(ip6)*scale(ip6))

      isL= INL(i-1) + 1
      ieL= INL(i  )
      !*voption indep (IAL,AL,SCALE)
      do k= isL, ieL
        inod= IAL(k)
        iq1= 6*inod - 5
        iq2= 6*inod - 5
        iq3= 6*inod - 3
        iq4= 6*inod - 2
        iq5= 6*inod - 1
        iq6= 6*inod

        AL(36*k-35)= AL(36*k-35)/(scale(ip1)*scale(iq1))
        AL(36*k-34)= AL(36*k-34)/(scale(ip1)*scale(iq2))
        AL(36*k-33)= AL(36*k-33)/(scale(ip1)*scale(iq3))
        AL(36*k-32)= AL(36*k-32)/(scale(ip1)*scale(iq4))
        AL(36*k-31)= AL(36*k-31)/(scale(ip1)*scale(iq5))
        AL(36*k-30)= AL(36*k-30)/(scale(ip1)*scale(iq6))

        AL(36*k-29)= AL(36*k-29)/(scale(ip2)*scale(iq1))
        AL(36*k-28)= AL(36*k-28)/(scale(ip2)*scale(iq2))
        AL(36*k-27)= AL(36*k-27)/(scale(ip2)*scale(iq3))
        AL(36*k-26)= AL(36*k-26)/(scale(ip2)*scale(iq4))
        AL(36*k-25)= AL(36*k-25)/(scale(ip2)*scale(iq5))
        AL(36*k-24)= AL(36*k-24)/(scale(ip2)*scale(iq6))

        AL(36*k-23)= AL(36*k-23)/(scale(ip3)*scale(iq1))
        AL(36*k-22)= AL(36*k-22)/(scale(ip3)*scale(iq2))
        AL(36*k-21)= AL(36*k-21)/(scale(ip3)*scale(iq3))
        AL(36*k-20)= AL(36*k-20)/(scale(ip3)*scale(iq4))
        AL(36*k-19)= AL(36*k-19)/(scale(ip3)*scale(iq5))
        AL(36*k-18)= AL(36*k-18)/(scale(ip3)*scale(iq6))

        AL(36*k-17)= AL(36*k-17)/(scale(ip4)*scale(iq1))
        AL(36*k-16)= AL(36*k-16)/(scale(ip4)*scale(iq2))
        AL(36*k-15)= AL(36*k-15)/(scale(ip4)*scale(iq3))
        AL(36*k-14)= AL(36*k-14)/(scale(ip4)*scale(iq4))
        AL(36*k-13)= AL(36*k-13)/(scale(ip4)*scale(iq5))
        AL(36*k-12)= AL(36*k-12)/(scale(ip4)*scale(iq6))

        AL(36*k-11)= AL(36*k-11)/(scale(ip5)*scale(iq1))
        AL(36*k-10)= AL(36*k-10)/(scale(ip5)*scale(iq2))
        AL(36*k-9 )= AL(36*k-9 )/(scale(ip5)*scale(iq3))
        AL(36*k-8 )= AL(36*k-8 )/(scale(ip5)*scale(iq4))
        AL(36*k-7 )= AL(36*k-7 )/(scale(ip5)*scale(iq5))
        AL(36*k-6 )= AL(36*k-6 )/(scale(ip5)*scale(iq6))

        AL(36*k-5 )= AL(36*k-5 )/(scale(ip6)*scale(iq1))
        AL(36*k-4 )= AL(36*k-4 )/(scale(ip6)*scale(iq2))
        AL(36*k-3 )= AL(36*k-3 )/(scale(ip6)*scale(iq3))
        AL(36*k-2 )= AL(36*k-2 )/(scale(ip6)*scale(iq4))
        AL(36*k-1 )= AL(36*k-1 )/(scale(ip6)*scale(iq5))
        AL(36*k   )= AL(36*k   )/(scale(ip6)*scale(iq6))
      enddo

      isU= INU(i-1) + 1
      ieU= INU(i  )
      !*voption indep (IAU,AU,SCALE)
      do k= isU, ieU
        inod= IAU(k)
        iq1= 6*inod - 5
        iq2= 6*inod - 4
        iq3= 6*inod - 3
        iq4= 6*inod - 2
        iq5= 6*inod - 1
        iq6= 6*inod

        AU(36*k-35)= AU(36*k-35)/(scale(ip1)*scale(iq1))
        AU(36*k-34)= AU(36*k-34)/(scale(ip1)*scale(iq2))
        AU(36*k-33)= AU(36*k-33)/(scale(ip1)*scale(iq3))
        AU(36*k-32)= AU(36*k-32)/(scale(ip1)*scale(iq4))
        AU(36*k-31)= AU(36*k-31)/(scale(ip1)*scale(iq5))
        AU(36*k-30)= AU(36*k-30)/(scale(ip1)*scale(iq6))

        AU(36*k-29)= AU(36*k-29)/(scale(ip2)*scale(iq1))
        AU(36*k-28)= AU(36*k-28)/(scale(ip2)*scale(iq2))
        AU(36*k-27)= AU(36*k-27)/(scale(ip2)*scale(iq3))
        AU(36*k-26)= AU(36*k-26)/(scale(ip2)*scale(iq4))
        AU(36*k-25)= AU(36*k-25)/(scale(ip2)*scale(iq5))
        AU(36*k-24)= AU(36*k-24)/(scale(ip2)*scale(iq6))

        AU(36*k-23)= AU(36*k-23)/(scale(ip3)*scale(iq1))
        AU(36*k-22)= AU(36*k-22)/(scale(ip3)*scale(iq2))
        AU(36*k-21)= AU(36*k-21)/(scale(ip3)*scale(iq3))
        AU(36*k-20)= AU(36*k-20)/(scale(ip3)*scale(iq4))
        AU(36*k-19)= AU(36*k-19)/(scale(ip3)*scale(iq5))
        AU(36*k-18)= AU(36*k-18)/(scale(ip3)*scale(iq6))

        AU(36*k-17)= AU(36*k-17)/(scale(ip4)*scale(iq1))
        AU(36*k-16)= AU(36*k-16)/(scale(ip4)*scale(iq2))
        AU(36*k-15)= AU(36*k-15)/(scale(ip4)*scale(iq3))
        AU(36*k-14)= AU(36*k-14)/(scale(ip4)*scale(iq4))
        AU(36*k-13)= AU(36*k-13)/(scale(ip4)*scale(iq5))
        AU(36*k-12)= AU(36*k-12)/(scale(ip4)*scale(iq6))

        AU(36*k-11)= AU(36*k-11)/(scale(ip5)*scale(iq1))
        AU(36*k-10)= AU(36*k-10)/(scale(ip5)*scale(iq2))
        AU(36*k-9 )= AU(36*k-9 )/(scale(ip5)*scale(iq3))
        AU(36*k-8 )= AU(36*k-8 )/(scale(ip5)*scale(iq4))
        AU(36*k-7 )= AU(36*k-7 )/(scale(ip5)*scale(iq5))
        AU(36*k-6 )= AU(36*k-6 )/(scale(ip5)*scale(iq6))

        AU(36*k-5 )= AU(36*k-5 )/(scale(ip6)*scale(iq1))
        AU(36*k-4 )= AU(36*k-4 )/(scale(ip6)*scale(iq2))
        AU(36*k-3 )= AU(36*k-3 )/(scale(ip6)*scale(iq3))
        AU(36*k-2 )= AU(36*k-2 )/(scale(ip6)*scale(iq4))
        AU(36*k-1 )= AU(36*k-1 )/(scale(ip6)*scale(iq5))
        AU(36*k   )= AU(36*k   )/(scale(ip6)*scale(iq6))
      enddo
    enddo

    deallocate(scale)
  end subroutine hecmw_solver_scaling_bk_66

end module hecmw_solver_scaling_66


text:
{#
  Copyright 2015 Andrew Dawson, Peter Dueben

  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at

      http://www.apache.org/licenses/LICENSE-2.0

  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.
#}
    ELEMENTAL FUNCTION {{ operator.name }}_{{ type1.name }}_{{ type2.name }} (x, y) RESULT (z)
        {{ type1.declaration }}, INTENT(IN) :: x
        {{ type2.declaration }}, INTENT(IN) :: y
        {{ operator.return_type.declaration }} :: z
        {% if operator.return_type.name == "rpe" %}
        z%sbits = MAX(significand_bits(x), significand_bits(y))
        {% endif %}
        z = x{{ type1.accessor }} {{ operator.operator }} y{{ type2.accessor }}
    END FUNCTION {{ operator.name }}_{{ type1.name }}_{{ type2.name }}

max_stars_repo_name: sparonuz/rpe

max_issues_repo_name: sparonuz/rpe

max_forks_repo_name: sparonuz/rpe


text:
	SUBROUTINE GSBARB  ( szbarb, ibrwid, ibrtyp, iret )
C************************************************************************
C* GSBARB								*
C* 									*
C* This subroutine sets the wind barb size, line width, and barb type.	*
C* If these parameters are not positive, no changes are made.		*
C*									*
C* NOTE THAT THE CALLING SEQUENCE FOR THIS SUBOUTINE HAS CHANGED FROM	*
C* PREVIOUS VERSIONS OF GEMPAK.						*
C* 									*
C* GSBARB  ( SZBARB, IBRWID, IBRTYP, IRET )				*
C*                                                                    	*
C* Input parameters:							*
C*	SZBARB		REAL		Barb size multiplier		*
C* 					  <=0 = no change		*
C*	IBRWID		INTEGER         Barb width multipier		*
C* 					  <=0 = no change		*
C*	IBRTYP		INTEGER		Barb type 			*
C*					  <=0 = no change		*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* I. Graffman/RDS	 5/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	 5/89	Make size = 0 no change			*
C* S. Schotz/GSC	 1/90   Added barb width and type parameters 	*
C* S. Schotz/GSC	 3/90   Fixed size to allow between 0 and 1	*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'DEVREQ.CMN'
	INCLUDE		'DEVSET.CMN'
C-----------------------------------------------------------------------
	iret = NORMAL
C
C*	First check if thes are the current requested characteristics.
C*	If so, do nothing.
C
	IF  ( szbarb .eq. rwbsz .and. ibrwid .eq. kbrwid .and.
     +        ibrtyp .eq. kbrtyp ) THEN
	ELSE
C
C*	    Set requested parameters
C
            IF  ( szbarb .gt. 0 ) rwbsz = szbarb
            IF  ( ibrwid .ge. 1 ) kbrwid = ibrwid
            IF  ( ibrtyp .ge. 1 ) kbrtyp = ibrtyp
C
C*	    Make changes only if the device has been set.
C
	    IF  ( ddev .ne. ' ' )  THEN
C
C*              Only set if requested different than what is already set
C
                IF ( rwbsz .ne. swbsz .or. kbrwid .ne. lbrwid .or.
     +               kbrtyp .ne. lbrtyp) THEN
		    CALL DSBARB  ( rwbsz, kbrwid, kbrtyp, swbsz,
     +                             lbrwid, lbrtyp, iret )
                END IF
	    END IF
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program era_sieve

implicit none
integer*4 :: N ! The natural number up to which prime numbers will be calculated
integer*4 :: ii, jj ! indexing variables
integer*4 :: allocate_status
logical, allocatable :: is_prime(:)
logical :: true = 1.gt.0, false = 0.gt.1

write(*,*) "Up to which number do you wish to get primes? Enter an integer N>=2."
read(*,*) N

! N should be at least 2
if (N.lt.2) stop "***The upper limit should be at least 2***"

allocate(is_prime(2:N),stat = allocate_status)
if (allocate_status .ne. 0) stop "***Not enough memory to allocate sieve vector***"
is_prime = (/(true, ii =2, N)/)

ii = 2

! Eratosthenes Sieve
do while (ii .le. sqrt(real(N,4))) ! Casting N as real to calculate sqrt(N)
	if (is_prime(ii)) then
		jj = ii**2
		do while (jj .le. N) ! Spotting multiples of ii starting from ii**2
			is_prime(jj) = false
			jj = jj + ii
		end do
	end if
	ii = ii + 1
end do

do ii = 2,N
	if (is_prime(ii)) then
		write(*,*) ii
	end if
end do

write(*,*) "Enter any number to normally exit the program."
read(*,*) N

end program era_sieve

max_stars_repo_name: eigen-carmona/quantum-computation-unipd

max_issues_repo_name: eigen-carmona/quantum-computation-unipd

max_forks_repo_name: eigen-carmona/quantum-computation-unipd


text:
ASUCD Senate Meeting ASUCD Senate/Minutes Minutes
October 19, 2006
5:30 Mee Room

I.  CALL TO ORDER

Vice President Genna Carnes called the meeting to order at 5:33 p.m.

II. QUORUM ROLL CALL

Present:  Carnes, Deepak, Farhad Motamed, Molnar, Rogers, RosasRomero, Salem, Sanders, Zamora, Markevich, Dardik, Khan, Xie, Lay, Rivera, Feldman

Absent:  Holloway, SavareeRuess, Ajlouny, Amaha, Higgins, Patel

III. STATUS OF PREVIOUS LEGISLATION PASSED

Sanders said all previous legislation has been signed.

IV. INTRODUCTION OF NEW LEGISLATION

Senate Resolution #2:  Authored by Harms, Introduced by Rogers, Referred to EAC
An ASUCD Senate Resolution to formally dissociate ASUCD from UCSA.

Senate Bill #14:  Authored by SavareeRuess, Introduced by Sanders, Referred to B&F
An ASUCD Senate Bill to amend the general program operating budget for 20062007.

Constitutional Amendment #19:  Authored by Ostrowski, Introduced by Ajlouny, Referred to
A Constitutional Amendment to change the necessary voting  5067

V. PUBLIC DISCUSSION

Sanders said Harney said that the majority of the Senate is defined in the constitution as a full senate, and all legislation needs to be passed by seven out of ten senators.  If not seven, even with abstinences, the legislation is not passed.  Also keep in mind that two thirds of the Senate equals nine, even if someone abstains.  This is the way it goes.  Also welcome to Sacramento New Technology High School, hanging out with us today.  At first when I saw all these people I was thinking “Oh! What did the senate do now?”

Hayley Steffen a former ECAC member and Alfredo Arredondo introduced selves as Student Assistants to the Chancellor.

Steffen said we want to make it a successful year and have an open agenda, we really want to collaborate with all of you.

Arredando said we have office hours every week, an hour before Senate meeting.  Also we will have cell phone numbers available on our business cards, so talk to us!

Zamora asked have you hired interns?

Steffen said we’re in the process of doing that.

Sanders moved to suspend to suspend the orders of the day and move into appointments and confirmations; properly seconded without objections.

VI. APPOINTMENTS AND CONFIRMATIONS

Sanders moved to confirm all members and/or alternates for EAC; properly seconded without objections.

Steven Lee, Amanda Peterson, Tiffany Setters, Dave Lue, Tai Gregor, Greg Dizon, and Luene Rizo introduced themselves.

Lay said they were all very impressive with good interviews.  It was a difficult decision, but the Senators and I all came to a consensus.  I urge a speedy confirmation.

Zamora asked why did you guys choose EAC?

Lee said in high school I did similar things working with the city of San Jose.

Peterson said similar reasons: I was part of a community department, I wish to establish ties between students and communities especially in Davis when community is such large part.  I also wish to raise student awareness.

Setters said I’m an out of state student and coming here has been extremely eye opening for me.  As corny as this may sound, I actually love Davis and have a hard time explaining my appreciation for the university as a whole.  I want to be a part of establishing strong ties between this university and the community.

Lue said I chose this commission because I understand that it works with legislation and voting issues.

Gregor said in high school I did a lot of lobbying.  Because of lobby corp and other things, there is a unique perception of this school around the state and even the country, and I think it’s important that it uphold its integrity.

Dizon said this commission looked most interesting.  I was involved in student government in high school and I just want to stay involved.

Rizo said sometimes when we’re in high school or college it’s easy for us to just hang out with our friends, worry about our own lives and we forget that there’s an entire world outside of school’s walls.  It’s important to remember that the community and everyone outside is valuable asset.

Deepak asked where are you all from?  What other type of outreaches or goals do you hope to accomplish (ie Keak da Sneak)?

Rizo said I’m from Roseville, California.  I used to work with the local government, both at the municipal and high school level (district office etc).  I think it’s important to be able to accomplish things from a government position; the UC is dependent upon state government: funding, legislation, etc. I also want to work with local business (ie charity concerts, record companies, entertainment).

Dizon said Huntington Beach (Orange County), California.  I’d like to see our EAC bringing in people other than entertainment for awareness in current issues (like AIDS, or really awareness of anything).  It’s nice to have concerts but also needed are programs to make students more aware of things affecting us today.

Gregor said I’m from Del Mar, California and my congressman is in jail.  It’s very important that our congressmen be in touch with people, and in touch with us as students.  It seems the state level is where congressman are more concerned about elderly issues (they have stronger voting) and we need to get on that.  Lobby corp and EAC wants representatives to address the issues that are important to our campus.  We need to visit them and make sure they’re listening to us.

Lue said Albany, California.  I want more people to turn out to vote.  Even my own friends don’t vote and I plan to hopefully turn that around. I think we just need to increase awareness of government on campus.

Setters said I’m from Eagle River, Alaska.  I really want stronger voter turnout, to be able to notify students about the benefits of voting, involvement, and funding positive relationships with community.  This positive attitude that’s already here really needs to continue to progress.

Peterson said Pollock Pines (near South Lake Tahoe).  I really want us to work on strengthening ties with local business and corporations.  Awareness extremely important; we can use programming to raise awareness.  An increased awareness and stronger relationship with the community has so many positive effects, including to increase voter registration and turnout.

Lee said in San Jose I was part of an organization called the Youth Advisory Council.  Through that I learned that on so many levels of government there are so many resources available but not utilized.  Some of these capital projects are geared toward our

demographic, so I hope all community/campus reach out to each other and be involved, give input.  When the community spends money it really should gear that toward our demographic.  Another issue is insurance and health at council level (we are the largest demographic uninsured).  State level school costs.  Help improve relationship with community versus the stereotype against college students.  Get goals accomplished!!

DO CONFIRM, 1200, EXTERNAL AFFAIRS COMMISSION: Steven Lee, Amanda Peterson, Tiffany Setters, Dave Lue, Tai Gregor, Greg Dizon, and Luene Rizo

Sanders moved to confirm any members and or alternates of EPPC; properly seconded without objection.

Emily mcael introduced himself.(check spelling)

Xie said she’s really passionate about this and gave a good interview.  I urge a speedy confirmation.

Sanders said I had the honor of sitting in on this interview. I remember you mentioning something about programs.  Why don’t you tell those to Senate?

Mcael said I wanted to join EPPC in order to be more aware of programs, and to help with making the campus more sustainable/renewable.  Also I hope that recycling programs can be further carried out, awareness of project compost needs to be increased and the program furthered.  Student awareness of environmental problems needs to be furthered and information given about how to live a more sustainable lifestyle.

DO CONFIRM, 1200, EMILY MCAEL, ENVIRONMENTAL POLICY AND PLANNING COMMISSION

Sanders moved to confirm any members and/or alternates of ECAC; properly seconded without objections.

Russell Manning introduced himself.

Khan said he was an excellent candidate and is experienced, he previously worked for the State of Black California’s government office.

Sanders said it sounds like you have great ideas what to do with campus, and the chair talked about your qualifications.  What direction would you like to take with the commission, and what ideas do you have for programs?

Manning said I want to create more communication between different ethnic groups, minority groups and ASUCD in general.  Definitely want to do anything that will create more forms of communication between groups.

DO CONFIRM, 1200, RUSSELL MANNING, ETHNIC AND CULTURAL AFFAIRS COMMISSION

Sanders motioned to resume the orders of the day; properly seconded without objections.

VII. PUBLIC ANNOUNCEMENTS

RosasRomero said it’s still Pride Week.  Tomorrow is Expressions From Within, a  celebration of personal art.  Staff/students will be bringing art to the LGBTRC, they’ll be on view from 1212:30 inside, plus free food, and people can just hang out.  From 12:302 we have a list of performances including one by Rob Roy!  Friday night at 8:00 in the  Veteran’s Memorial Building is the annual drag show: Davis is Burning!  Be there.

Xie said EPPC’s first event is in four minutes, MU 2, a public/urban sustainability discussion.  They’ll be discussing issues of growth/no growth in Davis, free food will be available, etc.

Markevich said UCSA has been a big talking point in AAC, many meetings discussing our relationship with them.  I’d like AAC to see Senate Resolution #2.

Vice President Carnes said this is public announcements, not discussion.  Bring it up later during discussion.

Feldman said Michelle Tea’s speaking engagement “The Personal is Still Political” is Thursday, October 26, 2006 at 1 p.m. in Voorhies Hall.

Vice President Carnes said CalPirg told me that officials are trying to take away early voting on campus; it’s something we usually have.  I only three contact info for city/county clerks.  If you care about this issue, shoot an email and make a call.  Don’t let them take away early voting.

Lay said next Wednesday, October 25th in the Coffee House will be a proposition debate, I’m contacting local student groups to speak on props 8589.  You should definitely be there to support us.

VIII. APPROVAL OF THE PAST MEETING’S MINUTES

Sanders said SGAO told me the minutes are available online, but they’re not there. Maybe they’ve been lost in cyber space?? Maybe Laabs knows where they are?

IX. UNIT DIRECTOR REPORTS
Vice President Carnes said we’ll come back to this when the unit director arrives.

X. EXOFICIO REPORTS

A.  Commission Chairs: Rivera, Lay, Markevich, Xie, Khan, Feldman, Dardik

B.  Executive: Carnes

XI. SENATE REPORTS

Sanders moved to suspend the orders of the day and return to unit director reports; properly seconded without objections.

XII. UNIT DIRECTOR REPORTS

Carlos Palacio (Aggie Student Store)

Patel asked what is this about candy being behind register?

Palacio said the store lost a considerable amount of money last year, this is my first quarter here so I’m left to guess why.  Seems that when candy is left in front of kids…

Molnar asked are you going…well…Sharon was talking about frozen dinners.  And what about microwaves and the Icee machines?

Palacio said today was the first day the Icee machine ran out of CO2.  We have burritos but can’t compete with the CoHo so we have to be selective.  Sharon at the CoHo is like my supervisor so I can’t compete with her.  Also we have a very small freezer.  But it is one of my main goals is to push the store in a direction to be more accessible for actual eating, not just snacking.

Amaha said you’re doing a good job.   I got the burrito today: it was good, organic too!

Palacio said it was good seeing you in the store.

Patel asked do you carry Fiji water?

Palacio said yes, 7up sells it.  It’s the most expensive water (because the bottle looks cool).

Markevich said it seems like the store gets overlooked between The Corral and Campus Copies.  How might you make it more visible?

Palacio said I had good idea: get a little cart to wheel around, put a sign on it and push around the Icee machine but maybe that would be made seasonal.  I was also thinking of putting seating and a table out front (ie bistro set), make store comfortable/loungy etc. so customers can chill out and talk to the…whatever.  We definitely need to do more to boost the sales to make this a profitable enterprise.

Lay asked do you sell It’s It icecream?

Palacio said we can look into it.

Molnar said are you going to sell baseball food (ie hotdogs)?

Palacio said at this point we’re starting simple, anything that I try to work into my system will have to require more employees.  I don’t have those resources to hire more pepple and especially with our tumultuous past I’m taking things one step at a time.  Caution is my priority. I want to familiarize myself with the functions of the store, be comfortable and be able to do everything with my eyes closed: at THAT point I will be ready to start implementing further dining systems.

Sanders moved to resume orders of day; properly seconded without objections.

XIII. OFFICER REPORTS

A. Senate: Sanders, Rogers, RosasRomero, Patel, Ajlouny, Deepak, Zamora,
    	     Amaha, Salem, Farhad Motamed, Molnar, Higgins

B. Media Board Rep:  Salem

C. Internet Networking Committee:  Zamora

XIV. NEW COURT CASES

Vice President Carnes said there are none.

XV. OTHER BUSINESS

Sanders said we have an interview coming up, need two senators.  It’s for the Whole Earth Festival Director Friday, October 20, 1:303.

Deepak said I’ll do it.

Patel said I can until 2, then RosasRomero will.

Vice President Carnes said to chairs: email Sanders if you want to schedule interviews, and give notice.  Don’t email whole table.

Sanders said it is my responsibility, but I’d like you to get it into me the Thursday before your interviews so I don’t have to call everyone.

XVI. PUBLIC DISCUSS

Markevich said AAC wants to see SR#2.

Vice President Carnes said it doesn’t apply to AAC, so I’d have to send it to everyone if I sent it to you.  We want you to keep in line with what commissions supposed to recommend, but I do appreciate your interest.

Ajlouny said sorry about my phone ringing, my brother is going into surgery and I’m waiting for mom to call me.

Higgins said I have a close friend who is throwing a Pleasure Party.  I’m not sure of the exact date/time, but come to me after the meeting and I will hit you up on facebook/call to let you know when and where.  It will be a lot of fun, you should check it out.  I’ll tell you what a pleasure party is afterwards.  Also George’s mom is the coolest.

RosasRomero said if you haven’t turned in the REACH application, they’re extended another week.  The director of the CCC told me that every senator was invited personally, there must’ve been a reason for that.

Vice President Carnes said I want to remind you to watch the language of this meeting, this is professional setting: especially considering that a high school could be visiting, etc.  Also if you didn’t make it to the training session or came late I have everything that was said transcribed down.  I’d be happy to have a discussion with you about what you missed.  We really should try to apply the things discussed, language, etc.

Sanders motioned to adjourn the meeting; properly seconded without objection.

XVII. ADJOURNMENT

Vice President Carnes adjourned the meeting at 6:40 p.m.

Submitted,

Bobby Gray
ASUCD Recorder


text:
Bus has multiple meanings.  You are probably looking for information on one of the following:

  Amtrak Bus.  The Amtrak passenger transportation system utilizes Amtrak Trains and Amtrak Buses.
  There no longer is a Greyhound Bus Station in Davis.  The Amtrak Buses stop at the Train Station.
  Davis Berkeley Shuttle is a bus providing daily weekday transport between the UC Davis and UC Berkeley campuses.
  Greyhound  A private, national bus line.
  Unitrans  Are the red buses that travel around Davis. All routes begin and end at UC Davis, either at the Memorial Union, the Silo or the Shields Library. Bus times are based on the UCD schedule and change during summer, winter and spring breaks and during finals.
  Yolobus  Yolobus is a bus system for travel throughout Yolo County and also to parts of Sacramento and Solano County Solano counties, including to the Sacramento International Airport.
  Electric bus  Sustainable transportation
  FAST  Fairfield and Suisun Transit

You may also want to check out Carfree in Davis and the UC Davis Commuter Challenge.

This is a Disambiguation disambiguation page  a navigational aid which lists other pages that might otherwise share the same title. If an article link referred you here, you might want to go back and fix it to point directly to the intended page.


text:
      integer e1,e4,e1p,e4p
      parameter(e1=5,e4=6,e1p=7,e4p=8)

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


    Program dgbbrd_example

!     DGBBRD Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_interfaces, Only: dgbbrd
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
      Character (1), Parameter :: vect = 'B'
!     .. Local Scalars ..
      Integer :: i, info, j, kl, ku, ldab, ldb, ldc, ldpt, ldq, m, n, ncc
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: ab(:, :), b(:, :), c(:, :), d(:), e(:), &
        pt(:, :), q(:, :), work(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, max, min
!     .. Executable Statements ..
      Write (nout, *) 'DGBBRD Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n, kl, ku, ncc
      ldab = kl + ku + 1
      ldb = m
      ldc = m
      ldpt = n
      ldq = m
      Allocate (ab(ldab,n), b(ldb,n), c(m,ncc), d(n), e(n-1), pt(ldpt,n), &
        q(ldq,m), work(2*m+2*n))

!     Read A from data file

      Read (nin, *)((ab(ku+1+i-j,j),j=max(i-kl,1),min(i+ku,n)), i=1, m)

!     Reduce A to upper bidiagonal form
      Call dgbbrd(vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, &
        c, ldc, work, info)

!     Print the absolute values of bidiagonal vectors d and e.
!     Any of these can differ by a sign change by combinations of sign
!     changes in columns of Q and P (rows of PT).
      Write (nout, *)
      Write (nout, *) 'Diagonal D:'
      Write (nout, 100) abs(d(1:n))
      Write (nout, *)
      Write (nout, *) 'Off-diagonal E:'
      Write (nout, 100) abs(e(1:n-1))
100   Format (1X, 4(3X,F11.4))

    End Program


text: #redirect Canada Geese


text:
program mpiTypes
use, intrinsic :: iso_fortran_env
use mpi_f08

implicit none

integer :: ierr, mrank, msize, vlen
character(MPI_MAX_LIBRARY_VERSION_STRING) :: version  ! allocatable not ok

print *,compiler_version()

call MPI_INIT(ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, mrank, ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, msize, ierr)
call MPI_GET_LIBRARY_VERSION(version, vlen, ierr)

call MPI_FINALIZE()

print '(A,I3,A,I3,A)', 'Image ', mrank, ' / ', msize, ':',version

print '(A25,A7)','type','value'
print '(A25,I4)','mpi_real',mpi_real
print '(A25,I4)','mpi_real8',mpi_real8

end program

max_stars_repo_name: cauachagas/sparse-fortran


!mm
      subroutine matrix_multiplication(A, B, n, C)

      implicit none

      real, intent(in) :: A(n, n)
      real, intent(in) :: B(n, n)
      real, intent(out) :: C(n, n)
      integer, intent(in) :: n

      integer :: i, j

      !$OMP PARALLEL DO
      do i = 1, n
          do j = 1, n
              C(j,i) = sum(A(:,j) * B(i,:))
          enddo
      enddo
      !$OMP END PARALLEL DO

      end subroutine matrix_multiplication


text:
      DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
      INTEGER          NEXT
      DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
      DATA   ZERO, ONE /0.0D0, 1.0D0/
C
C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
C     INCREMENT INCX .
C     IF    N .LE. 0 RETURN WITH RESULT = 0.
C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C           C.L.LAWSON, 1978 JAN 08
C
C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
C     HOPEFULLY APPLICABLE TO ALL MACHINES.
C         CUTLO = MAXIMUM OF  DSQRT(U/EPS)  OVER ALL KNOWN MACHINES.
C         CUTHI = MINIMUM OF  DSQRT(V)      OVER ALL KNOWN MACHINES.
C     WHERE
C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
C
C     BRIEF OUTLINE OF ALGORITHM..
C
C     PHASE 1    SCANS ZERO COMPONENTS.
C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C     VALUES FOR CUTLO AND CUTHI..
C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
C                   UNIVAC AND DEC AT 2**(-103)
C                   THUS CUTLO = 2**(-51) = 4.44089E-16
C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C                   THUS CUTHI = 2**(63.5) = 1.30438E19
C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
      DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C
      IF(N .GT. 0) GO TO 10
         DNRM2  = ZERO
         GO TO 300
C
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
      NN = N * INCX
C                                                 BEGIN MAIN LOOP
      I = 1
   20    GO TO NEXT,(30, 50, 70, 110)
   30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
      ASSIGN 50 TO NEXT
      XMAX = ZERO
C
C                        PHASE 1.  SUM IS ZERO
C
   50 IF( DX(I) .EQ. ZERO) GO TO 200
      IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
C
C                                PREPARE FOR PHASE 2.
      ASSIGN 70 TO NEXT
      GO TO 105
C
C                                PREPARE FOR PHASE 4.
C
  100 I = J
      ASSIGN 110 TO NEXT
      SUM = (SUM / DX(I)) / DX(I)
  105 XMAX = DABS(DX(I))
      GO TO 115
C
C                   PHASE 2.  SUM IS SMALL.
C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
   70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
C
C                     COMMON CODE FOR PHASES 2 AND 4.
C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
C
  110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
         SUM = ONE + SUM * (XMAX / DX(I))**2
         XMAX = DABS(DX(I))
         GO TO 200
C
  115 SUM = SUM + (DX(I)/XMAX)**2
      GO TO 200
C
C
C                  PREPARE FOR PHASE 3.
C
   75 SUM = (SUM * XMAX) * XMAX
C
C
C     FOR REAL OR D.P. SET HITEST = CUTHI/N
C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
C
   85 HITEST = CUTHI/FLOAT( N )
C
C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
C
      DO 95 J =I,NN,INCX
      IF(DABS(DX(J)) .GE. HITEST) GO TO 100
   95    SUM = SUM + DX(J)**2
      DNRM2 = DSQRT( SUM )
      GO TO 300
C
  200 CONTINUE
      I = I + INCX
      IF ( I .LE. NN ) GO TO 20
C
C              END OF MAIN LOOP.
C
C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
      DNRM2 = XMAX * DSQRT(SUM)
  300 CONTINUE
      RETURN
      END

ext: for


! 2017, S.A. Sato, Modified the format of the code
! ----------------------------------------------------------------------------
! Numerical diagonalization of 3x3 matrcies
! Copyright (C) 2006  Joachim Kopp
! ----------------------------------------------------------------------------
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Lesser General Public
! License as published by the Free Software Foundation; either
! version 2.1 of the License, or (at your option) any later version.
!
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with this library; if not, write to the Free Software
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
! ----------------------------------------------------------------------------

! ----------------------------------------------------------------------------
      SUBROUTINE ZHEEVC3(A, W)
! ----------------------------------------------------------------------------
! Calculates the eigenvalues of a hermitian 3x3 matrix A using Cardano's
! analytical algorithm.
! Only the diagonal and upper triangular parts of A are accessed. The access
! is read-only.
! ----------------------------------------------------------------------------
! Parameters:
!   A: The hermitian input matrix
!   W: Storage buffer for eigenvalues
! ----------------------------------------------------------------------------
! Dependencies:
!   SQRABS()
! ---------------------------------------------------------------------------
!     .. Arguments ..
      COMPLEX(8)       A(3,3)
      DOUBLE PRECISION W(3)

!     .. Parameters ..
      DOUBLE PRECISION SQRT3
      PARAMETER        ( SQRT3 = 1.73205080756887729352744634151D0 )

!     .. Local Variables ..
      DOUBLE PRECISION M, C1, C0
      COMPLEX(8)       DE
      DOUBLE PRECISION DD, EE, FF
      DOUBLE PRECISION P, SQRTP, Q, C, S, PHI

!     .. External Functions ..
      DOUBLE PRECISION SQRABS
      EXTERNAL         SQRABS

!     Determine coefficients of characteristic poynomial. We write
!           | A   D   F  |
!      A =  | D*  B   E  |
!           | F*  E*  C  |
      DE    = A(1,2) * A(2,3)
      DD    = SQRABS(A(1,2))
      EE    = SQRABS(A(2,3))
      FF    = SQRABS(A(1,3))
      M     = DREAL(A(1,1)) + DREAL(A(2,2)) + DREAL(A(3,3))
      C1    = ( DREAL(A(1,1)) * DREAL(A(2,2)) &
      + DREAL(A(1,1)) * DREAL(A(3,3)) &
      + DREAL(A(2,2)) * DREAL(A(3,3)) ) &
      - (DD + EE + FF)
      C0    = DREAL(A(3,3))*DD + DREAL(A(1,1))*EE + DREAL(A(2,2))*FF &
              - DREAL(A(1,1))*DREAL(A(2,2))*DREAL(A(3,3)) &
              - 2.0D0*(DREAL(A(1,3))*DREAL(DE)+DIMAG(A(1,3))*DIMAG(DE))

      P     = M**2 - 3.0D0 * C1
      Q     = M*(P - (3.0D0/2.0D0)*C1) - (27.0D0/2.0D0)*C0
      SQRTP = SQRT(ABS(P))

      PHI   = 27.0D0 * ( 0.25D0 * C1**2 * (P - C1) &
      + C0 * (Q + (27.0D0/4.0D0)*C0) )
      PHI   = (1.0D0/3.0D0) * ATAN2(SQRT(ABS(PHI)), Q)

      C     = SQRTP * COS(PHI)
      S     = (1.0D0/SQRT3) * SQRTP * SIN(PHI)

      W(2) = (1.0D0/3.0D0) * (M - C)
      W(3) = W(2) + S
      W(1) = W(2) + C
      W(2) = W(2) - S

      END SUBROUTINE
! End of subroutine ZHEEVC3


text:
      subroutine setCh(n)
      include 'dswacom.h'

      integer n,k
      do  k=1, 29
        wabr(k)=0
        wabr(n)=1
      enddo
      return
      end

      subroutine printCh
      include 'dswacom.h'
      integer n,k
      if(wabr(1).gt.0) write(*,*) "  S1 S1            ",wabr(1)
      if(wabr(2).gt.0) write(*,*) "  S1 S2            ",wabr(2)
      if(wabr(3).gt.0) write(*,*) "  S2 S2            ",wabr(3)
      if(wabr(4).gt.0) write(*,*) "  S3 S3            ",wabr(4)
      if(wabr(5).gt.0) write(*,*) "  S1 S3            ",wabr(5)
      if(wabr(6).gt.0) write(*,*) "  S2 S3            ",wabr(6)
      if(wabr(7).gt.0) write(*,*) "  S- S+            ",wabr(7)
      if(wabr(8).gt.0) write(*,*) "  Z S1             ",wabr(8)
      if(wabr(9).gt.0) write(*,*) "  Z S2             ",wabr(9)
      if(wabr(10).gt.0) write(*,*) "  Z S3             ",wabr(10)
      if(wabr(11).gt.0) write(*,*) "  W- S+ and W+ S-  ",wabr(11)
      if(wabr(12).gt.0) write(*,*) "  Z0 Z0            ",wabr(12)
      if(wabr(13).gt.0) write(*,*) "  W+ W-            ",wabr(13)
      if(wabr(14).gt.0) write(*,*) "  nu_e nu_e-bar    ",wabr(14)
      if(wabr(15).gt.0) write(*,*) "  e+ e-            ",wabr(15)
      if(wabr(16).gt.0) write(*,*) "  nu_mu nu_mu-bar  ",wabr(16)
      if(wabr(17).gt.0) write(*,*) "  mu+ mu-          ",wabr(17)
      if(wabr(18).gt.0) write(*,*) "  nu_tau nu_tau-bar",wabr(18)
      if(wabr(19).gt.0) write(*,*) "  tau+ tau-        ",wabr(19)
      if(wabr(20).gt.0) write(*,*) "  u u-bar          ",wabr(20)
      if(wabr(21).gt.0) write(*,*) "  d d-bar          ",wabr(21)
      if(wabr(22).gt.0) write(*,*) "  c c-bar          ",wabr(22)
      if(wabr(23).gt.0) write(*,*) "  s s-bar          ",wabr(23)
      if(wabr(24).gt.0) write(*,*) "  t t-bar          ",wabr(24)
      if(wabr(25).gt.0) write(*,*) "  b b-bar          ",wabr(25)
      if(wabr(26).gt.0) write(*,*) "  gluon gluon      ",wabr(26)
      if(wabr(27).gt.0) write(*,*) "  q q gluon (not im",wabr(27)
      if(wabr(28).gt.0) write(*,*) "  gamma gamma (1-lo",wabr(28)
      if(wabr(29).gt.0) write(*,*) "  Z0 gamma (1-loop)",wabr(29)
      end

      subroutine setMcdmSdSI(Mcdm,csSD,csSI)
      include 'dswacom.h'
      real*4 Mcdm,csSI,cdSD
       write(*,*) "test",  Mcdm,csSD,csSI

       wamwimp=Mcdm
       wasigsip=csSI
       wasigsdp=csSD
      end

      subroutine printMcdmSdSI
      include 'dswacom.h'
      write(*,*) '!!!', wamwimp,wasigsip,wasigsdp
      end


text:
Brain Maps is a UC Davis Locally Owned Websites website with extremely high resolution scans of brains.

Edward Jones is the Principal Investigator for the project.