text:
subroutine ThetaTauLeap(K,T,Y0,Y,options,info)
    !
    !   Copyright 2016, Viktor Reshiak, All rights reserved.
    !
    !   Purpose
    !   =======
    !   Find solution of a_tau_s stochastic chemical system.
    !
    !
    !   Method
    !   ======
    !   Theta Tau-leaping:
    !                                                          /                                                   \
    !  / y1'[k] \    / y1[k] \    / nu_11  nu_12 ... nu_1m \   |        / prop_1'[k] \               / prop_1[k] \ |
    !  | y2'[k] |    | y2[k] |    | nu_21  nu_22 ... nu_2m |   |        | prop_2'[k] |               | prop_2[k] | |
    !  |    .   |  = |   .   |  + |   .     .          .   | * | eta1 * |     .      |  + (1-eta1) * |   .       | | * (1-theta) * tau
    !  |    .   |    |   .   |    |   .           .    .   |   |        |     .      |               |   .       | |
    !  \ yn'[k] /    \ yn[k] /    \ nu_n1  nu_n2 ... nu_nm /   |        \ prop_n'[k] /               \ prop_n[k] / |
    !                                                          \                                                   /
    !
    !                                                          /                                 \
    !  / y1''[k] \   / y1'[k] \   / nu_11  nu_12 ... nu_1m \   | / P_1' \   / prop_1'[k] \       |
    !  | y2''[k] |   | y2'[k] |   | nu_21  nu_22 ... nu_2m |   | | P_2' |   | prop_2'[k] |       |
    !  |    .    | = |   .    | + |   .     .          .   | * | |  .   | - |     .      | * tau |
    !  |    .    |   |   .    |   |   .           .    .   |   | |  .   |   |     .      |       |
    !  \ yn''[k] /   \ yn'[k] /   \ nu_n1  nu_n2 ... nu_nm /   | \ P_m' /   \ prop_n'[k] /       |
    !                                                          \                                 /
    !   where P_j' = Poisson( prop_j(y'[k]) * tau )
    !                                                           /                                                      \
    !  / y1[k+1] \   / y1''[k] \   / nu_11  nu_12 ... nu_1m \   |        / prop_1[k+1] \               / prop_1''[k] \ |
    !  | y2[k+1] |   | y2''[k] |   | nu_21  nu_22 ... nu_2m |   |        | prop_2[k+1] |               | prop_2''[k] | |
    !  |   .     | = |   .     | + |   .     .          .   | * | eta2 * |     .       |  + (1-eta2) * |     .       | | * theta * tau
    !  |   .     |   |   .     |   |   .           .    .   |   |        |     .       |               |     .       | |
    !  \ yn[k+1] /   \ yn''[k] /   \ nu_n1  nu_n2 ... nu_nm /   |        \ prop_n[k+1] /               \ prop_n''[k] / |
    !                                                           \                                                      /
    !
    !
    !   IN
    !   ==
    !   1) num_species         - dimension of the state vector
    !   2) num_reactions         - number of reactions
    !   3) K         - number of time points
    !   4) prop      - function to propensity vector, its Jacobian and rates of reactions subsystems
    !   5) nu        - num_species-by-num_reactions stoichiometric matrix
    !   6) T         - K-dimensional vector of time points
    !   7) Y0        - num_species-dimensional vector with initial data
    !   8) options   - (optional) struct with options of the solver
    !       theta    - implicitness parameter
    !       Jacobian - true if Jacobain is supplied with 'prop' function, false otherwise
    !
    !
    !   OUT
    !   ===
    !   Y            - num_species-by-K solution array. Each row in Y is the solution of the corresponding equation
    !   info         - structure with information about solvers
    !     time       - average time of solver per time step
    !     iter       - average number of nonlinear iterations per time step
    !
    !
    !   Notes
    !   =====
    !   It is assumed that time grid is uniform

    implicit none

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! dummy arguments

    integer(kind=4), intent(in)                               :: K
       real(kind=8), intent(in),  dimension(num_species)      :: Y0
       real(kind=8), intent(in),  dimension(K)                :: T
       real(kind=8), intent(out), dimension(:,:), allocatable :: Y

    type(solver_options), intent(in),  optional    :: options
    type(solver_info),    intent(out), optional    :: info

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    type(solver_options)         :: op
    type(solver_info)            :: inf
    type(fsolve_info)            :: fsolve_info

    real(kind=8)                 :: tau, theta
    real(kind=8), dimension(num_reactions) :: a_tau_s
    integer,      dimension(num_reactions) :: dN
    real(kind=8), dimension(num_reactions,num_species) :: Ja
    real(kind=8), dimension(num_species,num_species) :: eye_N

    integer(kind=8)  :: count1, count2, count_rate
    integer(kind=4)  :: i, skipped
    integer(kind=4)  :: my_thread

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if ( present(options) ) then
        op = options
    else
        op = default_solver_options
    endif

    if ( .not. allocated(Y) ) then
        allocate(Y(num_species,K))
    endif

    tau   = T(2) - T(1)
    theta = op%theta

    eye_N = 0.d0
    forall(i=1:num_species) eye_N(i,i) = 1.d0

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    my_thread = omp_get_thread_num()+1

    inf%iter      = 0
    inf%neg_steps = 0
    inf%residual  = 0.d0
    call system_clock(count1)
    Y(:,1) = Y0
	!do i = 2,K
    i = 2
    skipped = 0
    do while (i <= K)

        ! stochastic step
        a_tau_s = propensities(Y(:,i-1))
        a_tau_s = a_tau_s * tau
        call poissrnd( num_reactions, a_tau_s, dN )
        Y(:,i) = Y(:,i-1) + matmul( nu, dN - a_tau_s * theta )

        call deterministic_corrector()

        ! update solution (preserve integer states)
        Y(:,i) = Y(:,i-1) + matmul( nu, nint( dN + ( propensities(Y(:,i))*tau - a_tau_s ) * theta ) );

        call preserve_positivity()

        i = i + 1
        if ( inf%iter == default_fsolve_options%MaxIter ) then
            i = i - 1;
            skipped = skipped + 1
            if ( skipped > K * 0.5 ) then
                print*, 'Solution did not converge!', dN, Y(:,i)
                ! need this for script to keep the execution window active
                print*, 'Press any key'
                read(*,*)
                call exit(0)
            endif
            print*, 'Skip iteration', inf%iter, skipped, i, K
        endif
	enddo
    call system_clock(count2, count_rate)
    inf%time     = (count2-count1) / real(count_rate,8)
    inf%residual = inf%residual / (K-1);
    !inf%iter     = inf%iter     / (K-1);

    if ( present(info) ) then
        info = inf
    endif

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    contains

    subroutine deterministic_corrector()
        theta_tau(:,my_thread) = theta * tau
        Y_buf(:,my_thread)     = Y(:,i)

        call fsolve( num_species, num_species, fun, Y(:,i), X0=Y(:,i-1), info=fsolve_info, Jac=Jac )

        inf%iter     = inf%iter     + fsolve_info%iter
        inf%residual = inf%residual + fsolve_info%final_residual
    end subroutine deterministic_corrector

    subroutine preserve_positivity()
        integer(kind=4) :: j, n

        if ( any( Y(:,i) < 0 ) ) then
            inf%neg_steps = inf%neg_steps + 1
            do j = 1,num_species
                do while ( Y(j,i)<0 )
                    do n = 1,num_reactions
                        if ( nu(j,n)<0 ) then
                            Y(:,i) = Y(:,i) - nu(:,n)
                        endif
                    enddo
                enddo
            enddo
        endif
    end subroutine preserve_positivity

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! objective function for deterministic predictor & corrector
    subroutine fun(dim_f,dim_X,X,f)
        integer,      intent(in)                     :: dim_f, dim_X
        real(kind=8), intent(in),  dimension(dim_X ) :: X
        real(kind=8), intent(out), dimension(dim_f)  :: f
        integer :: my_thread

        my_thread = omp_get_thread_num()+1

        f = X - Y_buf(:,my_thread) - matmul( nu, propensities(X) * theta_tau(:,my_thread) )
    end subroutine fun

    ! objective Jacobian for deterministic predictor & corrector
    subroutine Jac(dim_f,dim_X,X,J)
        integer,      intent(in)                          :: dim_f, dim_X
        real(kind=8), intent(in),  dimension(dim_X)       :: X
        real(kind=8), intent(out), dimension(dim_f,dim_X) :: J
        integer :: my_thread

        my_thread = omp_get_thread_num()+1

        Ja = propJacobian(X)

        forall(i=1:num_species) Ja(:,i) = Ja(:,i) * theta_tau(:,my_thread)

        J = eye_N - matmul( nu, Ja )
    end subroutine Jac

end subroutine ThetaTauLeap

max_stars_repo_name: revitmt/ChemicalKinetics

max_issues_repo_name: revitmt/ChemicalKinetics

max_forks_repo_name: revitmt/ChemicalKinetics


text:
C23456789012345678901234567890123456789012345678901234567890123456789012
C
c      Program to compute spline fits to fermi integrals
cc  Must provide data file 94
      subroutine initferm(FilePath)

!     use parallel_module
!     use mpi

      implicit  double precision (a-h,o-z)
      parameter (n=201)
      character*128 FilePath
      dimension f32(n),f12(n),fm12(n),eta(n),fr(n)
      dimension f32a(n),f12a(n),fra(n),fia(n)
      common /spl/eta,f32,f12,fr,f32a,f12a,fra,fia
      dimension send_buf(4,n)
C
!     if (myid.eq.0) then
      open(94,file=TRIM(FilePath)//'/fermi.atb',
     & status='old')
C
      do 10 i=1,n
       read(94,*)eta(i),f32(i),f12(i),fm12(i)
       send_buf(1,i) = eta(i)
       send_buf(2,i) = f32(i)
       send_buf(3,i) = f12(i)
       send_buf(4,i) = fm12(i)
 10   continue
C
      close(94,status='keep')

!     endif

!     i_extent=n*4
!     CALL MPI_BCAST( send_buf , i_extent, MPI_DOUBLE_PRECISION, 0,
!    &                     MPI_COMM_WORLD, ierr)

      do 20 i=1,n
       eta(i) = send_buf(1,i)
       f32(i) = send_buf(2,i)
       f12(i) = send_buf(3,i)
       fm12(i) = send_buf(4,i)
 20    fr(i)=f12(i)/fm12(i)

C
       call spline(eta,f12,n,f12a)
       call spline(eta,f32,n,f32a)
       call spline(eta,fr,n,fra)
       call spline(f12,eta,n,fia)
       return
       end subroutine initferm

ext: f

max_stars_repo_name: srichers/weaklib


text:
      SUBROUTINE DN2LRD(DR, IV, L, LH, LIV, LV, ND, NN, P, R, RD, V)
C
C  ***  COMPUTE REGRESSION DIAGNOSTIC AND DEFAULT COVARIANCE MATRIX FOR
C        DRN2G  ***
C
C  ***  PARAMETERS  ***
C
      INTEGER LH, LIV, LV, ND, NN, P
      INTEGER IV(LIV)
      DOUBLE PRECISION DR(ND,P), L(LH), R(NN), RD(NN), V(LV)
C
C  ***  CODED BY DAVID M. GAY (WINTER 1982, FALL 1983)  ***
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      DOUBLE PRECISION DD7TPR
      EXTERNAL DD7TPR, DL7ITV, DL7IVM,DO7PRD, DV7SCP
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER COV, I, J, M, STEP1
      DOUBLE PRECISION A, FF, S, T
C
C  ***  CONSTANTS  ***
C
      DOUBLE PRECISION NEGONE, ONE, ONEV(1), ZERO
C
C  ***  INTRINSIC FUNCTIONS  ***
C/+
      DOUBLE PRECISION DSQRT
C/
C
C  ***  IV AND V SUBSCRIPTS  ***
C
      INTEGER F, H, MODE, RDREQ, STEP
C/6
C     DATA F/10/, H/56/, MODE/35/, RDREQ/57/, STEP/40/
C/7
      PARAMETER (F=10, H=56, MODE=35, RDREQ=57, STEP=40)
C/
C/6
C     DATA NEGONE/-1.D+0/, ONE/1.D+0/, ZERO/0.D+0/
C/7
      PARAMETER (NEGONE=-1.D+0, ONE=1.D+0, ZERO=0.D+0)
C/
      DATA ONEV(1)/1.D+0/
C
C++++++++++++++++++++++++++++++++  BODY  +++++++++++++++++++++++++++++++
C
      STEP1 = IV(STEP)
      I = IV(RDREQ)
      IF (I .LE. 0) GO TO 999
      IF (MOD(I,4) .LT. 2) GO TO 30
      FF = ONE
      IF (V(F) .NE. ZERO) FF = ONE / DSQRT(DABS(V(F)))
      CALL DV7SCP(NN, RD, NEGONE)
      DO 20 I = 1, NN
         A = R(I)**2
         M = STEP1
         DO 10 J = 1, P
            V(M) = DR(I,J)
            M = M + 1
 10         CONTINUE
         CALL DL7IVM(P, V(STEP1), L, V(STEP1))
         S = DD7TPR(P, V(STEP1), V(STEP1))
         T = ONE - S
         IF (T .LE. ZERO) GO TO 20
         A = A * S / T
         RD(I) = DSQRT(A) * FF
 20      CONTINUE
C
 30   IF (IV(MODE) - P .LT. 2) GO TO 999
C
C  ***  COMPUTE DEFAULT COVARIANCE MATRIX  ***
C
      COV = IABS(IV(H))
      DO 50 I = 1, NN
         M = STEP1
         DO 40 J = 1, P
            V(M) = DR(I,J)
            M = M + 1
 40         CONTINUE
         CALL DL7IVM(P, V(STEP1), L, V(STEP1))
         CALL DL7ITV(P, V(STEP1), L, V(STEP1))
         CALL DO7PRD(1, LH, P, V(COV), ONEV, V(STEP1), V(STEP1))
 50      CONTINUE
C
 999  RETURN
C  ***  LAST LINE OF DN2LRD FOLLOWS  ***
      END

max_stars_repo_name: DanIverson/OpenVnmrJ

max_issues_repo_name: timburrow/openvnmrj-source

max_forks_repo_name: timburrow/openvnmrj-source


text:
HANSARD REVISE * NUMERO 3
Le mercredi 24 septembre 1997
L'AVION D'ENTRAINEMENT A REACTION MONITOR
LE DISCOURS DU TRONE
M. Guy Saint-Julien
L'UNIVERSITE DU QUEBEC A CHICOUTIMI
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
L'hon. Jean J. Charest
Le tres hon. Jean Chretien
L'hon. Jean J. Charest
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
L'hon. Pierre S. Pettigrew
Le tres hon. Jean Chretien
Le tres hon. Jean Chretien
LE DISCOURS DU TRONE
L'hon. Pierre S. Pettigrew
L'hon. Pierre S. Pettigrew
LE DIRECTEUR GENERAL DES ELECTIONS DU CANADA
La Loi sur la sante et la loi sur les drogues
LE DISCOURS DU TRONE
Reprise du debat sur l'Adresse
Le tres hon. Jean Chretien
L'hon. Jean J. Charest
HANSARD REVISE * NUMERO 3
Le mercredi 24 septembre 1997
La seance est ouverte a 14 heures.
Tous les deputes de levent et chantent o Canada.]
L'AVION D'ENTRAINEMENT A REACTION MONITOR
M. Bob Wood (Nipissing, Lib.):
J'espere que tous les deputes seront la ce soir.
M. Inky Mark (Dauphin-Swan River, Ref.):
Je ne serais pas ici aujourd'hui si cette loi n'avait pas ete abrogee.
L'annee 1997 marque le 50e anniversaire de l'abrogation de cette loi.
Mme Judi Longfield (Whitby-Ajax, Lib.):
Mme Pierrette Venne (Saint-Bruno-Saint-Hubert, BQ):
M. Stan Keyes (Hamilton-Ouest, Lib.):
Mme Elinor Caplan (Thornhill, Lib.):
Le Canada prend enfin les moyens necessaires.
Comme cette question est particulierement importante pour la communaute juive...
Le depute de Nanaimo-Cowichan.
M. Reed Elley (Nanaimo-Cowichan, Ref.):
C'est avec beaucoup de tristesse que les Canadiens ont appris son deces.
M. Gary Pillitteri (Niagara Falls, Lib.):
Cet important evenement a renforce la cooperation qui existe entre nos deux grands pays.
M. Daniel Turp (Beauharnois-Salaberry, BQ):
LE DISCOURS DU TRONE
Mme Marlene Jennings (Notre-Dame-de-Grace-Lachine, Lib.):
Je vous accorde encore quelques secondes.
M. Guy Saint-Julien (Abitibi, Lib.):
Monsieur le President, la sante de l'economie canadienne ne fait plus de doute.
Voici de bonnes nouvelles pour les Canadiens et les Canadiennes.
Mme Angela Vautour (Beausejour-Petitcodiac, NPD):
Le gouvernement ne respecte pas sa propre Loi sur l'equite salariale.
M. Keith Martin (Esquimalt-Juan de Fuca, Ref.):
Nous avons tous pleure.
Elle nous a montre par l'exemple comment surmonter nos difficultes.
Diana a defendu la cause des enfants malades, des sideens et des grands malades incurables.
L'UNIVERSITE DU QUEBEC A CHICOUTIMI
M. Denis Coderre (Bourassa, Lib.):
Mme Elsie Wayne (Saint John, PC):
Les Canadiens ages sont en droit de compter sur une retraite assuree.
Les jeunes Canadiens veulent que le RPC leur soit accessible lorsqu'ils en auront besoin.
Les Canadiens a revenu moyen seront les plus durement touches.
Mme Albina Guarnieri (Mississauga-Est, Lib.):
Aucun etre humain n'a autant accompli pour autant de gens, et pour si peu.
On se souviendra d'elle simplement et affectueusement sous le nom de Mere.
Ainsi, les declarations preliminaires de ceux qui posent des questions seront plus breves.
Je demanderais egalement aux ministres de collaborer en donnant des reponses plus courtes.
Alors, nous allons voir ce qui se passera aujourd'hui.
Nous surveillerons la duree des interventions.
M. Preston Manning (chef de l'opposition, Ref.):
S'agissait-il la d'une simple omission?
Le tres hon. Jean Chretien (premier ministre, Lib.):
Le budget est prevu pour le mois de fevrier.
M. Preston Manning (chef de l'opposition, Ref.):
Le tres hon. Jean Chretien (premier ministre, Lib.):
M. Preston Manning (chef de l'opposition, Ref.):
Le tres hon. Jean Chretien (premier ministre, Lib.):
Il essaie de trouver une cible sur laquelle il puisse tirer.
M. Monte Solberg (Medicine Hat, Ref.):
Pourtant, selon son ministere, c'est de la foutaise.
Qui a raison et qui a tort?
L'hon. Paul Martin (ministre des Finances, Lib.):
Monsieur le President, c'est le depute qui a tort.
M. Monte Solberg (Medicine Hat, Ref.):
Les Canadiens veulent une reduction d'impot.
Je voudrais une reponse de l'aspirant premier ministre, s'il vous plait.
L'hon. Paul Martin (ministre des Finances, Lib.):
Monsieur le President, il est bon que le depute fasse partie de l'opposition officielle.
Pour la premiere fois, il comprend les realites politiques du pays.
Permettez-moi de donner quelques chiffres.
Sous le gouvernement actuel...
M. Gilles Duceppe (Laurier-Sainte-Marie, BQ):
Pourtant, la realite est tout autre.
Le tres hon. Jean Chretien (premier ministre, Lib.):
Cela prouve que nous pouvons travailler en partenariat a l'interieur du Canada.
M. Gilles Duceppe (Laurier-Sainte-Marie, BQ):
On considere les provinces comme de grosses municipalites.
Le tres hon. Jean Chretien (premier ministre, Lib.):
Prenons, par exemple, les revendications traditionnelles du Quebec sur la main-d'oeuvre.
Nous avons conclu une entente avec les gouvernements provinciaux, y inclus celui du Quebec.
M. Michel Gauthier (Roberval, BQ):
Ma question s'adresse au premier ministre.
Le tres hon. Jean Chretien (premier ministre, Lib.):
Et je suis tres content de l'attitude des premiers ministres a ce sujet.
M. Michel Gauthier (Roberval, BQ):
Le tres hon. Jean Chretien (premier ministre, Lib.):
Mme Alexa McDonough (Halifax, NPD):
Monsieur le President, j'ai attendu ce moment longtemps.
Ma question s'adresse au premier ministre.
Le gouvernement l'a fait en ce qui concerne la reduction du deficit.
Le tres hon. Jean Chretien (premier ministre, Lib.):
Je suis heureux qu'elle soit parmi nous.
Nous devons le faire pour creer des emplois.
C'est de cette facon que nous allons creer des emplois.
Mme Alexa McDonough (Halifax, NPD):
Le tres hon. Jean Chretien (premier ministre, Lib.):
Permettez-moi de rappeler certains commentaires que j'ai lus le 10 fevrier dernier.
L'hon. Jean J. Charest (Sherbrooke, PC):
Il pourrait poser un geste des aujourd'hui en reduisant les cotisations d'assurance-emploi.
Le tres hon. Jean Chretien (premier ministre, Lib.):
L'hon. Jean J. Charest (Sherbrooke, PC):
Pour plus de clarte, je pose au premier ministre une question directe et simple.
Le tres hon. Jean Chretien (premier ministre, Lib.):
Il n'avait jamais eu d'excedent dans le fonds d'assurance-chomage.
A ce chapitre, nous sommes au premier rang des pays du G7.
Mme Deborah Grey (Edmonton-Nord, Ref.):
Le tres hon. Jean Chretien (premier ministre, Lib.):
Voila exactement ce que nous entendons par la.
Mme Deborah Grey (Edmonton-Nord, Ref.):
Monsieur le President, le concept de l'egalite et le mot egalite sont importants.
Le tres hon. Jean Chretien (premier ministre, Lib.):
C'est tres important.
L'egalite veut qu'on reconnaisse les gens pour ce qu'ils sont.
Malgre l'egalite, nous sommes aux prises avec un enorme probleme de pauvrete au Canada.
Il faut donc que les riches aident les pauvres.
Le depute de Temiscamingue.
M. Pierre Brien (Temiscamingue, BQ):
Monsieur le President, ma question s'adresse au premier ministre.
M. Pierre Brien (Temiscamingue, BQ):
L'hon. Pierre S. Pettigrew (ministre du Developpement des ressources humaines, Lib.):
Il s'agit d'une enquete longitudinale nationale que nous faisons deja depuis quelques annees.
M. Rahim Jaffer (Edmonton-Strathcona, Ref.):
Le tres hon. Jean Chretien (premier ministre, Lib.):
Monsieur le President, je remercie le depute pour sa question.
M. Rahim Jaffer (Edmonton-Strathcona, Ref.):
Le tres hon. Jean Chretien (premier ministre, Lib.):
Il y a eu un referendum.
Le choix des senateurs par voie electorale etait au nombre des reformes que nous preconisions.
Le Parti reformiste a vote contre cet accord.
Pour l'instant, les provinces se penchent sur la question.
Le moment venu, nous agirons.
Nous ne voulons pas refaire ce qui se fait deja dans les autres provinces.
LE DISCOURS DU TRONE
M. Yvan Loubier (Saint-Hyacinthe-Bagot, BQ):
Ma question s'adresse au ministre des Finances.
L'hon. Paul Martin (ministre des Finances, Lib.):
M. Yvan Loubier (Saint-Hyacinthe-Bagot, BQ):
L'hon. Paul Martin (ministre des Finances, Lib.):
M. Grant Hill (Macleod, Ref.):
L'hon. Allan Rock (ministre de la Sante, Lib.):
Il faudra pour cela travailler tres fort.
Le depute de Macleod.
M. Grant Hill (Macleod, Ref.):
Ils n'en ont pas la moindre idee.
L'assurance-maladie a ete touchee et est mal en point.
L'hon. Allan Rock (ministre de la Sante, Lib.):
M. Paul Crete (Kamouraska-Riviere-du-Loup-Temiscouata-Les Basques, BQ):
Monsieur le President, ma question s'adresse au ministre du Developpement des ressources humaines.
L'hon. Pierre S. Pettigrew (ministre du Developpement des ressources humaines, Lib.):
Notamment, nous avons annonce une strategie de l'emploi-jeunesse qui est extremement importante.
Nous travaillons tres fort du cote du chomage, et vous le savez.
Cette strategie, qui cree des experiences de travail, fonctionne bien.
M. Paul Crete (Kamouraska-Riviere-du-Loup-Temiscouata-Les Basques, BQ):
L'hon. Pierre S. Pettigrew (ministre du Developpement des ressources humaines, Lib.):
Nous avons mis sur pied un fonds transitoire de creation d'emplois.
Nous ne croyons pas au support passif du revenu.
C'est ce que veulent nos travailleurs.
M. Carmen Provenzano (Sault Ste. Marie, Lib.):
Monsieur le President, ma question s'adresse au secretaire parlementaire du ministre des Affaires etrangeres.
M. Ted McWhinney (secretaire parlementaire du ministre des Affaires etrangeres, Lib.):
Elles sont probablement contraires aux dispositions de l'ALENA.
M. John Duncan (Ile de Vancouver-Nord, Ref.):
En fait, le gouvernement gere mal ce dossier depuis quatre ans.
L'hon. David Anderson (ministre des Peches et des Oceans, Lib.):
Monsieur le President, l'affirmation du depute est inexacte.
La question du traite du saumon est certes importante.
M. John Duncan (Ile de Vancouver-Nord, Ref.):
Le gouvernement continue de negliger la Colombie-Britannique.
L'hon. David Anderson (ministre des Peches et des Oceans, Lib.):
LE DIRECTEUR GENERAL DES ELECTIONS DU CANADA
LA LOI SUR LA SANTE ET LA LOI SUR LES DROGUES
M. Jim Gouk (West Kootenay-Okanagan, Ref.):
M. John Finlay (Oxford, Lib.):
Monsieur le President, je vous felicite pour votre reelection.
J'ai deux petitions a presenter aujourd'hui.
M. John Finlay (Oxford, Lib.):
Monsieur le President, la deuxieme petition est signee par 57 electeurs de ma circonscription.
M. Alex Shepherd (Durham, Lib.):
Monsieur le President, je demande que toutes les questions demeurent au Feuilleton .
LE DISCOURS DU TRONE
REPRISE DU DEBAT SUR L'ADRESSE
M. Preston Manning (chef de l'opposition, Ref.):
Je tiens a vous feliciter pour votre reelection a la presidence.
Au nom des deputes de l'opposition officielle, je vous souhaite la meilleure des chances.
Felicitations, monsieur, au nom de nous tous.
Nous felicitons egalement la deputee de Halifax et nos collegues neo-democrates.
Je felicite le chef du Parti conservateur et ses collegues.
Je felicite aussi le premier ministre.
Il a une longue carriere en politique.
Nous avons bien tente d'y mettre fin, mais ca n'a pas marche.
Nous aurions prefere etre a sa place.
Il n'appartient pas a mon parti.
Il appartient a mes electeurs.
A mon avis, leur message est quadruple.
C'est signe de mecontentement.
Il y a un message la-dedans.
Il nous incombe de nous assurer qu'il en soit ainsi.
Le troisieme message est des plus interessant.
Le Bloc est exclusivement un parti quebecois.
Nous considerons avoir un mandat a deux volets.
Je veux feliciter le gouvernement lorsqu'il le merite.
Nous voulons critiquer les aspects de ce discours que nous jugeons inacceptables.
Ainsi, les gens peuvent se demander dans quelle voie nous allons nous engager.
Dans le discours du Trone, on pretend regarder en avant.
Il est important de voir la realite qui se cache derriere les mots.
Il etait a bout de souffle.
Macdonald, leur guide, n'etait plus la.
Le discours du trone cite Sir Wilfrid Laurier.
Il parlait de la situation financiere.
Le Parti liberal, qui fut autrefois un grand parti, est en train de s'epuiser.
Nous voulons examiner les os secs.
Nous voulons rendre a Cesar ce qui appartient a Cesar.
Examinons tout d'abord l'aspect economique.
Le gouvernement propose bien peu a cet egard.
Il ne mentionne aucun secteur de ressources naturelles.
Il ne dit rien au sujet du secteur manufacturier.
Il ne traite pas de l'investissement dans la connaissance et la creativite.
La question est de savoir quoi faire avec les trois autres pneus.
Ou sont les principes qui nous guideront dans ces domaines?
Permettez-moi de toucher un mot sur la securite sociale.
C'est aussi notre avis.
Le filet de securite sociale exige, selon nous, une nouvelle approche.
Elle devait faire quelque chose pour resserrer la Loi sur les jeunes contrevenants.
Elle devait revoir le processus de determination de la peine.
Elle devait modifier le systeme de liberation conditionnelle.
Ou est la reaction aux besoins qu'ont exprimes les victimes de crimes?
Ces gens ont toute notre sympathie.
Pour terminer, j'aborderai maintenant la question inepuisable de l'unite nationale.
C'etait la position du gouvernement en 1995.
Je le repete, nous assistons maintenant a quelques changements modestes.
Ces changements sont beaucoup trop modestes.
C'est bien maigre.
Permettez-moi de donner quelques exemples de ces propositions.
Le passage du discours du Trone consacre a l'unite est extremement bref.
Nous n'avons rien contre.
Pourquoi prendre une telle chance quand il est question d'unite nationale?
Il pourrait s'y prendre de deux facons.
Cela a entache les negociations de 1982.
Ensuite, il y a eu les negociations de Charlottetown.
La population pense qu'il s'est agi d'une consultation apres coup.
En quoi a consiste notre experience?
Pendant un certain nombre d'annees, nous avons laisse M. Trudeau conduire la voiture.
En retrospective, nous lui avons pardonne.
Que leur a repondu M. Trudeau?
Nous avons ensuite eu un autre conducteur, le premier ministre Mulroney.
Il a dit que, pour cela, il fallait se rendre a Charlottetown.
Nous y sommes alles et nous ne sommes toujours pas parvenus a un Canada uni.
Un peu plus tard, notre actuel premier ministre prend le volant.
Il conduit tres tres lentement.
Assis a l'arriere, nous lui demandons si nous sommes arrives a un Canada uni.
Nous avons trouve de quoi remplir deux pages.
Nous promulguerions des lois pour renforcer l'union economique canadienne.
Ou est la mesure legislative visant a faciliter l'etablissement de normes nationales?
Il n'y a rien dans ce texte qui touche vraiment aux problemes des autochtones.
Et ou est le projet de loi parallele?
Cela donnerait a la Chambre la meilleure representation regionale qu'elle ait jamais eue.
Qu'en est-il de la question de la representation proportionnelle dans cette enceinte?
Ou est le projet de loi qui traiterait de cette question?
La Chambre est censee prevoir.
Tout cela doit se faire a l'avance.
Il faut d'abord adopter une loi sur les assemblees constituantes.
Le Canada disposerait alors des outils pour faciliter les choses.
Nous aurions une nouvelle repartition des pouvoirs pour le XXIe siecle.
Il n'y a rien de tout cela dans le discours du Trone.
Difficile d'imaginer meilleur moyen d'entraver la reprise economique.
Et qu'en est-il du Nord, nos dernieres terres vierges, demographiquement et ecologiquement parlant?
Il n'y a aucune reconnaissance du Nord, aucune vision pour le Nord.
Meme Diefenbaker avait une vision pour le Nord.
J'ai garde le meilleur pour la fin: l'Ouest.
Je prierais les deputes ministeriels d'ecouter un instant.
Le parti reformiste est le principal porte-parole de l'Ouest nouveau au Parlement federal.
L'Ouest peut trouver lui-meme la solution a ses griefs regionaux.
Il en fait d'ailleurs la preuve tous les jours.
Cet exemple peut etre reproduit en maints autres endroits au Canada.
Je termine en commentant une fois de plus un discours du Trone sterile.
Peut-etre certains se contenteraient de pas grand-chose.
Nous devons parler au nom de ceux que le gouvernement n'a pas ecoutes.
Nous devons proposer des politiques que le gouvernement n'a pas elaborees.
Pour resumer, je propose:
Le tres hon. Jean Chretien (premier ministre, Lib.):
Je leur predis un tres bel avenir a la Chambre des communes.
J'ai ete elu depute pour la onzieme fois.
Ils m'ont appris que les gens sont ce qui compte le plus en politique.
Nos agriculteurs font leurs recoltes et les jeunes reprennent le chemin de l'ecole.
Nous avons renoue avec notre vitalite economique.
Nous investirons dans les enfants, notre ressource la plus precieuse.
Notre pays a investi dans l'assurance-maladie il y a exactement 30 ans.
Je siegeais alors a la Chambre.
Le tres hon. Jean Chretien:
Desormais, les problemes actuels ne doivent plus etre legues aux generations futures.
Une nouvelle vague d'optimisme deferle sur le Canada.
Nous devons maintenant progresser ensemble vers un nouveau millenaire.
Beaucoup de deputes sont des nouveaux venus a la Chambre.
Certains se plaisent a denigrer le Parlement, mais ils ont tort.
La situation actuelle est bien meilleure que lorsque je suis devenu depute.
Pour y parvenir, j'ai travaille avec des collegues des deux cotes de la Chambre.
Nous avons tous partage ce succes.
Au cours de ce mandat, notre gouvernement assurera la securite des aines pour l'avenir.
Permettez-moi de dire en toute franchise que nous avons beaucoup de travail a faire.
C'est, tout au contraire, un gouvernement capable de mieux s'en acquitter.
La population a egalement besoin de notre aide.
Le chomage chez les jeunes est trop eleve.
Meilleure sera l'instruction des jeunes, meilleures seront leurs chances de trouver un emploi.
Et c'est precisement ce que nous allons faire.
Le savoir et l'apprentissage comptent parmi les plus importants de ces domaines.
La technologie est un univers de promesses, mais elle doit avoir une ame.
Trop souvent, nous oublions ou nous ignorons tout ce que nous avons accompli ensemble.
En travaillant ensemble, Baldwin et Lafontaine nous auront donne un gouvernement responsable.
Nous avons bati un pays etonnant, respecte et meme envie du monde entier.
Le ministre des Affaires etrangeres merite toutes nos felicitations.
Ce traite sauvera des milliers de vies humaines.
Les Canadiens ont manifeste cette force lors des inondations au Saguenay et au Manitoba.
Notre gouvernement ressent cette force morale dans tout ce qu'il entreprend.
C'est un message qui doit etre entendu.
C'est un autre pas important pour batir la comprehension mutuelle et la confiance.
Il ne faudrait pas y voir plus.
Et notre ennemi, c'est la confusion.
L'Alberta et la Saskatchewan n'etaient pas encore des provinces.
Pourtant, nous savions que l'avenir s'ouvrait devant nous.
Au debut du siecle, Laurier exprimait ce reve d'avenir lorsqu'il a declare:
La cathedrale etait faite de granit, de chene et de marbre.
Nous avons bati ce pays et nous continuons d'en faconner les elements.
Les jeunes le feront a leur tour au prochain siecle.
Leur architecture sera nouvelle mais demeurera canadienne.
Cette grandeur aura acquis un sens nouveau mais demeurera canadienne.
Notre pays sera un modele pour le monde.
Nous avons tous le tres grand honneur de sieger au Parlement canadien.
Pour bien des pays, le Canada est le pays reve.
M. Gilles Duceppe (Laurier-Sainte-Marie, BQ):
Voila ou ce deuxieme constat, par son objectif, rejoint notre premier constat.
Cela veut dire aussi plus de professeurs, plus de services pour les etudiants.
Ensuite, il faut utiliser les surplus pour mener la lutte contre la pauvrete.
Le pire empietement jamais vu dans l'histoire du Canada.
Meme Pierre Elliott Trudeau n'a jamais ose se rendre aussi loin.
C'est donc peu dire.
A quand les normes nationales dans l'education?
A quand les examens nationaux dans l'education?
Une attitude qui se reflete de facon manifeste dans le discours du Trone du gouvernement.
Manifestement, le scenario de Charlottetown risque fort de se reproduire.
C'est ca, le caractere bilingue du Canada.
Le Bloc quebecois n'a jamais deroge de ses principes de base.
Plus que jamais, nous jugeons que la souverainete est necessaire.
L'impasse constitutionnelle persiste toujours et encore plus depuis la declaration de Calgary.
On assiste a une veritable surenchere sur le dos du peuple quebecois.
Voila la cause de tout le blocage.
Les Canadiens forment un peuple.
Les autochtones forment des peuples.
Les gens partout a travers le monde forment des peuples.
En conclusion, permettez-moi quelques reflexions sur l'avenir du Canada.
Il nous invite meme a feter par anticipation cette reussite a venir.
Notre lecture de la realite est tout autre.
Selon nous, le Canada est a la croisee des chemins.
C'est la voie de la raison, c'est la voie du bon sens.
Cette voie mene au cul-de-sac.
En terminant, je propose:
Qu'on modifie l'amendement en ajoutant, apres les mots programme legislatif ce qui suit:
L'amendement est recevable.
Nous disposons maintenant de quelques minutes pour des questions et commentaires.
Mme Alexa McDonough (Halifax, NPD):
C'est sans aucun doute un honneur, mais egalement une tache intimidante.
Je veux aussi feliciter et remercier les autres candidats a la Presidence.
C'est pour nous une celebration a la fois douce et amere.
Ces livres presentaient une formule plus democratique pour le Canada.
On dit souvent de Halifax que c'est un peu une circonscription barometre.
Mais, a ma connaissance, personne ne recommande pareille formule d'equilibre des sexes.
Cela aere l'esprit, calme les nerfs et purge l'ame.
Aggraver nos inegalites ne servira qu'a rejouir ceux qui cherchent a nous diviser.
Ce Parlement doit etre axe sur la reconstruction de ce Canada.
Les femmes sont les plus vulnerables.
Il n'est pas surprenant que le gouvernement ait abandonne les femmes.
Les familles et les enfants du Canada attendent encore.
Voyons comment le gouvernement traite ses propres fonctionnaires feminins.
Que se passe-t-il vraiment avec l'assurance-maladie?
Qui en fait les frais?
Ou va le gouvernement dans le dossier environnemental?
Il convient toutefois d'etre justes.
Il y a un engagement gouvernemental ferme dont nous nous rejouissons.
Je veux rendre a Cesar ce qui revient a Cesar.
Tous les Canadiens peuvent etre fiers de cette realisation.
J'imagine que le chef du Nouveau Parti democratique est en train de recapituler.
Monsieur le President, je voudrais demander le consentement unanime pour conclure.
Je sais que le temps est limite et je respecte cela.
L'hon. Jean J. Charest (Sherbrooke, PC):
Monsieur le President, cela me rappelle la premiere journee de mon retour a la Chambre.
Je me disais que j'etais entoure de progressistes conservateurs.
Quel contraste avec la derniere legislature.
Je dois vous dire que je suis tres tres endette envers ces gens.
J'ai eu le grand privilege de les representer depuis plusieurs annees.
C'est precisement ce que la derniere campagne nous a permis de faire.
Cela a toujours ete le cas, d'ailleurs.
Mais ils ont une memoire bien selective.
Notre pays affiche egalement une excellente performance au chapitre du commerce international.
Dieu merci, nous avons conclu des accords commerciaux.
Nous nous sommes battus et avons conclu l'Accord de libre-echange nord-americain.
Or, que voyons-nous aujourd'hui?
Il est maintenant tout fier de parler de ses nouveaux  amigos, au Chili.
Il veut conclure des accords avec l'Argentine.
Il veut en conclure avec Israel.
Il n'en a jamais assez.
Et pourtant, ils nous ont livre une lutte acharnee sur cette question.
Notre economie interieure est mal en point depuis quatre ans.
Qu'applaudissaient-ils encore?
Quoi d'autre encore?
Peut-etre devraient-ils aller voir ce qui se passe dans le Canada atlantique.
Une reduction sans precedent de 40 p. 100 du financement.
Y a-t-il eu la moindre consultation?
Le systeme canadien d'assurance maladie est mal en point.
Les liberaux devront en repondre au cours de la presente legislature.
Nous nous attendions a un programme pour le pays.
Qu'avons-nous eu a la place?
Les voeux pieux habituels.
Combien de fois faut-il redecouvrir le probleme du chomage chez les jeunes?
Environ un demi-million de jeunes Canadiens sont au chomage.
Je crois que leur nombre depasse les 400 000.
Quelle grande initiative du gouvernement liberal.
Plusieurs messages ont ete lances durant la campagne.
Le premier ministre et moi sommes d'accord sur au moins l'un d'eux.
Il y a donc un message qui leur a ete envoye.
Qu'a-t-on eu a la place?
Les deputes voudront le remarquer, dans quels programmes le gouvernement reinvestira-t-il?
On reinvestira dans la culture, apres avoir coupe.
On reinvestira dans l'education postsecondaire, qu'ils ont coupee.
Il y a quelques minutes, on a vu un spectacle extraordinaire.
Cela se passait en 1995.
Et savez-vous ce qui est arrive, monsieur le President?
Cet apres-midi, le premier ministre annonce un nouveau programme de bourses.
Je les reconnais bien.
Mais qu'est-ce qui a change entre-temps?
Il n'est pas trop tard.
Permettez-moi de dire un secret au premier ministre.
Sa position sur l'assurance maladie est insoutenable.
Cela ne fonctionnera pas.
Nous sommes impatients de voir ce partenariat dont parle le premier ministre.
Le gouvernement national peut faire quelque chose d'utile a cet egard.
Vous n'avez pas besoin de demander la permission a qui que ce soit.
Tout ce que vous avez a faire, c'est travailler avec les provinces...
Je sais que le depute a beaucoup d'experience a la Chambre.
L'hon. Jean J. Charest:
Je serai tres heureux de m'adresser a la presidence, monsieur le President.
Comme je le disais, le premier ministre a l'occasion d'agir.
Nous esperons qu'il est sincere lorsqu'il se dit pret a accepter les conseils.
La partie la plus decevante a trait aux emplois et aux impots.
Le gouvernement aurait pu faire quelque chose.
Le gouvernement fonctionne en suivant quelques principes constants.
C'est la philosophie liberale.
C'est plutot bien tourne.
C'est comme cela que les choses se passent.
Le gouvernement aurait pu organiser les choses differemment pour le pays.
Je ne suis peut-etre pas convaincant.
Peut-etre suis-je trop sectaire.
Les deputes pensent-ils que la Chambre de commerce du Canada est trop sectaire?
Pensent-ils que le Conseil du patronat est trop sectaire?
Je ne le crois pas.
J'espere que le gouvernement finira par apprendre.
Y a-t-il unanimite pour que le depute termine ses observations?
L'hon. Jean J. Charest:
Monsieur le President, le gouvernement parle d'impots.
Il parle de commerce.
Il ne cesse de parler de retablir l'equilibre de la federation.
On va reequilibrer les responsabilites dans la federation.
Encore une fois, je donnerai un conseil a mon collegue.
Le premier ministre effectuera un autre voyage avec Equipe Canada en janvier.
Puis-je lui faire une humble suggestion?
Qui s'y opposerait?
C'est aussi une des raisons qui justifient fortement ce leadership.
Voila donc certaines des idees auxquelles nous croyons.
La reforme des pensions sera aussi tres importante.
Le premier ministre sait que les Canadiens sont tres inquiets et ils devraient l'etre.
Le gouvernement a un objectif secret.
La classe moyenne va se faire avoir par le gouvernement.
Les Canadiens de la classe moyenne seront touches.
Nous voulons que le gouvernement dise les choses comme elles sont.
Nous sommes prets a debattre de la reforme, mais nous voulons la verite.
C'est une charge sociale.
On ne peut les remplacer.
Et on reagit a ce qu'ils font.
Personne ne s'attend non plus a ce qu'on livre des amendements constitutionnels.
Ils n'ont pas a se poser la question, ils le sont, point.
On va avancer dans la mesure ou on peut avancer ensemble.
C'en est une des lecons.
Mais notre premier ministre, malheureusement, ne semble pas enclin a le faire.
C'est ce sur quoi on jugera ce Parlement.
L'honorable whip du Parti reformiste invoque le Reglement.
Je pense qu'il est temps de conclure.
Nous avons deja depasse la periode prescrite de 15 a 20 minutes.
Je vous prierais de bien vouloir lui demander de conclure.
Le depute n'a pas invoque le Reglement.
Il ne fait que presenter son point de vue.
L'hon. Jean J. Charest:
Monsieur le President, le depute pourrait partir a sa guise.
Il n'est pas tenu de rester assis et de supporter cela.
Je ne veux pas m'eterniser la-dessus.
Si ca ne lui plait pas, il n'a qu'a quitter les lieux.
Voila une chose sur laquelle nous, le Parlement, avons un effet des plus direct.
(La seance est levee a 18 h 17.)


text:
PROGRAM hash_test
  use hash_chtch_class

  character(len=128)::res

  type(hash_chtch)   :: h1,h2

  h1 = hash_chtch()

  print*, '----'
  CALL h1%hash_set("one",  '1.0')
  CALL h1%hash_set("two",  '2.0')
  CALL h1%hash_set("three",'3.0')

  print*, '----'
  CALL h1%hash_print
  CALL h1%hash_get("one", res)
  PRINT*, res
  CALL h1%hash_get("two", res)
  PRINT*, res
  CALL h1%hash_get("three", res)
  PRINT*, res

  !print*, '----'
  !CALL h1%hash_set("one", '60.0')
  !CALL h1%hash_print
  !CALL h1%hash_get("one", res)
  !PRINT*, res
END PROGRAM hash_test

max_stars_repo_name: LukeMcCulloch/FortranHashTables

max_issues_repo_name: LukeMcCulloch/FortranHashTables

max_forks_repo_name: LukeMcCulloch/FortranHashTables


! { dg-do run }

module reduction5
  intrinsic ior, min, max
end module reduction5

  call test1
  call test2
contains
  subroutine test1
    use reduction5, bitwise_or => ior
    integer :: n
    n = Z'f'
!$omp parallel sections num_threads (3) reduction (bitwise_or: n)
    n = ior (n, Z'20')
!$omp section
    n = bitwise_or (Z'410', n)
!$omp section
    n = bitwise_or (n, Z'2000')
!$omp end parallel sections
    if (n .ne. Z'243f') call abort
  end subroutine
  subroutine test2
    use reduction5, min => max, max => min
    integer :: m, n
    m = 8
    n = 4
!$omp parallel sections num_threads (3) reduction (min: n) &
!$omp & reduction (max: m)
    if (m .gt. 13) m = 13
    if (n .lt. 11) n = 11
!$omp section
    if (m .gt. 5) m = 5
    if (n .lt. 15) n = 15
!$omp section
    if (m .gt. 3) m = 3
    if (n .lt. -1) n = -1
!$omp end parallel sections
    if (m .ne. 3 .or. n .ne. 15) call abort
  end subroutine test2
end

! { dg-final { cleanup-modules "reduction5" } }


text:
! -*- F90 -*-

subroutine commoninit()
  include 'common.inc'
  integer i

  if (commoninitflag .ne. 'commonblockinitdone') then
     !print *, "Initialising LHAPDF steering data"
     commoninitflag = 'commonblockinitdone'

     ! LHAPDF common block
     lhaname = ' '
     lhaset = 0
     lhamemb = 0

     ! LHASETS common block
     do i = 1, nmxset
        lhanames(i) = ' '
        lhanumbers(i) = 0
        lhamembers(i) = 0
     end do
     nsets = 0

     ! LHAPDFC common block
     lhapath = 'pdfsets'

     ! LHACONTROL common block

     do i = 1, 20
        lhaparm(i) = ' '
        lhavalue(i) = 0.0d0
     end do

     ! LHAGLSTA common block
     xminnum = 0.0d0
     xmaxnum = 0.0d0
     q2minnum = 0.0d0
     q2maxnum = 0.0d0
     totnum = 0.0d0
     xminnup = 0.0d0
     xmaxnup = 0.0d0
     q2minnup = 0.0d0
     q2maxnup = 0.0d0
     totnup = 0.0d0
  end if

end subroutine commoninit

max_stars_repo_name: AllaMaevskaya/AliRoot

max_issues_repo_name: AllaMaevskaya/AliRoot

max_forks_repo_name: AllaMaevskaya/AliRoot


text:
      subroutine Root_Mover()
      Include 'public.ins'
      Include 'puplant.ins'
      Double precision A,B,C,P,Sum
      Character InString*132
      Integer newjjj
C  AD added RMass, RMassM, MRL into Public.ins 2011-08-01
C AD 12-2-2011 this code from MaizSim10.11 is based on Yakov Pachepsky's diffusion code

      Dimension A(MBandD,NumNPD),B(NumNPD),F(NumNPD),DS(NumNPD),
     &            Gc(NumNPD),Sc(NumNPD),Fc(NumNPD),
     &            MRL(NumNPD),RDenTotal(NumNPD)
      Common /RootM/  DMolx,DMolz,Vel,Ac(NumNPD),
     !                NLevel,dt,epsi,CourMax,lUpW,
     !                S(3,3),Wz(3),Wx(3),Ri(3),Ci(3),Bi(3),List(3),
c     !         VxROld(NumNPD),VzROld(NumNPD),
c     !         VxR(NumNPD), VzR(NumNPD), hOld(NumNPD),tOld(NumNPD),
c     !         VxRH(NumNPD),VzRH(NumNPD),WeTab(3,2*NumElD),
     !                Dispzz(NumNPD),Dispxx(NumNPD),Dispxz(NumNPD),
     !                RMassOld(NumNPD),RDenTotal
      If (lInput.eq.0) goto 11

      Open(40,file = VarietyFile, status = 'old',ERR = 10)

      im=70
      il=0
55       Read (40,'(A132)') InString
          if (InString(1:14).ne.'[RootDiff]') goto 55
      im=im+1
      il=il+1
      Read(40,*,ERR=10)
      im=im+1
      il=il+1
      Read(40,*,ERR=10)
      im=im+1
      il=il+1
      Read(40,*,ERR=10) epsi,lUpW,CourMax
      Nlevel=1
      if(epsi.lt.0.999) NLevel=2
      im=im+1
      il=il+1
      Read(40,*,ERR=10)
      im=im+1
      il=il+1
      Read(40,*,ERR=10) DMolx,DMolz,Vel
      close(40)

      Do n=1,NumNP
        RMassOld(n)=RMassY(n)
c        VxR(n)=0.0D0
c        VxROld(n)=0.0D0
c        VzROld(n)=VzR(n)
      Enddo
      Movers(4)=1
      Return
C
C  Routine calculations
11    Continue

C  No Plant = No Root Activity
      If(NShoot.eq.0) Return

      t=Time
      dt = Step
      xMul=1.0D0
      alf=1.0D0-epsi
      newjjj=MBand
c
      Do 13 n=1,NumNP
       DS(n)=0.0D0
       Gc(n)=0.0D0
       Sc(n)=0.0D0
       B(n) =0.0D0
       Fc(n)=0.0D0
c          Ac(i)=0.
         If(lOrt) B1(n)=RMassY(n)

          If(epsi.lt.0.001D0) then
            if(lOrt) newjjj = IADD(n)
            A(newjjj,n)=0.0D0
          Else

          Do 12 j=1,MBandD
              A(j,n)=0.0D0
12          Continue
          Endif
13      Continue
c
c        A loop to calculate geotropic velocity
CDT     we are going to not use this as a first approximation.
c
c        Do n=1,NumNP
c           VzR(n)=Vel*f1(hNew(n))*f2(Tmpr(n))
c     &         *f3(RDenY(n)+RDenM(n))
C          VzR(i)=Vel*RGCF(i)
C        Enddo
C
C  assign fluxes of root growth
        Do n=1,NumNP
          Gc(n)=ADWR(n)/nodeArea(n)
        Enddo
C
C Assembling matrixes
C
        Do 21 Level=1,NLevel
         If(lUpW.eq.1) then
!          If(Level.eq.NLevel) then
!            Do n=1,NumNP
!cccz test there is no convection
!                if ((abs(VxR(n))+abs(VzR(n))).gt.0.0D0) then
!                 aaaa=1
!                endif
!cccz
!             VxRH(n)=VxR(n)
!             VzRH(n)=VzR(n)
!            Enddo
!          Else
!	      Do n=1,NumNP
!cccz test there is no convection
!             if ((abs(VxROld(n))+abs(VzROld(n))).gt.0.0D0) then
!                 aaaa=1
!             endif
!cccz
!             VxRH(n)=VxROld(n)
!             VzRH(n)=VzROld(n)
!            Enddo
!          Endif

          Call RDisper(NumNP,Dispxx,Dispzz,Dispxz,DMolx,DMolz,
     &             hNew,Tmpr,RDenTotal)
!          Call RWeFact(NumNP,NumEl,NumElD,x,y,KX,WeTab,VxRH,VzRH,
!     &             Dispxx,Dispzz,Dispxz)
         Endif
         If(Level.eq.NLevel) then
          do 14 n=1,NumNP
           Ac(n)=-1
14        continue
         Endif
         Do 15 n=1,NumNP
            F(n)=0.
            If(Level.eq.NLevel) DS(n)=0.
15        Continue
          NumSEl=0
          Do 19 n=1,NumEl
            NUS=4
            If(KX(n,3).eq.KX(n,4)) NUS=3
            Do 18 k=1,NUS-2
              NumSEl=NumSEl+1
              i=KX(n,1)
              j=KX(n,k+1)
              l=KX(n,k+2)
              List(1)=i
              List(2)=j
              List(3)=l
              Ci(1)=x(l)-x(j)                    ! cccz the order: ijl is counter-clockwise
              Ci(2)=x(i)-x(l)
              Ci(3)=x(j)-x(i)
              Bi(1)=y(j)-y(l)
              Bi(2)=y(l)-y(i)
              Bi(3)=y(i)-y(j)
              AE=(Ci(3)*Bi(2)-Ci(2)*Bi(3))/2.
              If(KAT.eq.1) xMul=2.*3.1416*(x(i)+x(j)+x(l))/3.
              FMul=xMul*AE/4.
              GcE=(Gc(i)+Gc(j)+Gc(l))/3.
              Ec1=(Dispxx(i)+Dispxx(j)+Dispxx(l))/3.
              Ec2=(Dispxz(i)+Dispxz(j)+Dispxz(l))/3.
              Ec3=(Dispzz(i)+Dispzz(j)+Dispzz(l))/3.
              AcE=(Ac(i)+Ac(j)+Ac(l))/3.
!              FcE=(Fc(i)+Fc(j)+Fc(l))/3.
!              VxE=(VxR(i)+VxR(j)+VxR(l))/3.
!              VzE=(VzR(i)+VzR(j)+VzR(l))/3.
              SMul1=-1./AE/4.*xMul
!              SMul2=AE/20.*xMul
!cccz test there is no convection
!               if ((abs(VxE)+abs(VzE)).gt.0.0D0) then
!                 aaaa=1
!               endif
!cccz
!              If(lUpW.eq.1) then
!                NS=NumSEl
!                W1=WeTab(1,NS)
!                W2=WeTab(2,NS)
!                W3=WeTab(3,NS)
!                Wx(1)=
!     &          2.*VxRH(i)*(W2-W3)+VxRH(j)*(W2-2.*W3)+VxRH(l)*(2.*W2-W3)
!                Wx(2)=
!     &          VxRH(i)*(2.*W3-W1)+2.*VxRH(j)*(W3-W1)+VxRH(l)*(W3-2.*W1)
!                Wx(3)=
!     &          VxRH(i)*(W1-2.*W2)+VxRH(j)*(2.*W1-W2)+2.*VxRH(l)*(W1-W2)
!                Wz(1)=
!     &          2.*VzRH(i)*(W2-W3)+VzRH(j)*(W2-2.*W3)+VzRH(l)*(2.*W2-W3)
!                Wz(2)=
!     &          VzRH(i)*(2.*W3-W1)+2.*VzRH(j)*(W3-W1)+VzRH(l)*(W3-2.*W1)
!                Wz(3)=
!     &          VzRH(i)*(W1-2.*W2)+VzRH(j)*(2.*W1-W2)+2.*VzRH(l)*(W1-W2)
!cccz cccz test there is no convection
!             if ((abs(VxRH(i))+abs(VxRH(j))+abs(VxRH(k))).gt.0.0D0) then
!                 aaaa=1
!              endif
!             if ((abs(VzRH(i))+abs(VzRH(j))+abs(VzRH(k))).gt.0.0D0) then
!                 aaaa=1
!         endif
!              Endif
            do 17 j1=1,3
              i1=List(j1)
              F(i1)=F(i1)+FMul*(GcE+Gc(i1)/3.)
              if(Level.eq.NLevel) DS(i1)=DS(i1)+FMul*(AcE+Ac(i1)/3.)
              do 16 j2=1,3
                i2=List(j2)
                S(j1,j2)=SMul1*(Ec1*Bi(j1)*Bi(j2)+Ec3*Ci(j1)*Ci(j2)+
     !                         Ec2*(Bi(j1)*Ci(j2)+Ci(j1)*Bi(j2)))
!                S(j1,j2)=S(j1,j2)-(Bi(j2)/8.*(VxE+VxR(i1)/3.)+
!     !                            Ci(j2)/8.*(VzE+VzR(i1)/3.))*xMul
!                if(lUpW.eq.1) S(j1,j2)=S(j1,j2)-xMul*
!     !                            (Bi(j2)/40.*Wx(j1)+Ci(j2)/40.*Wz(j1))
!                ic=1
!                if(i1.eq.i2) ic=2
!                S(j1,j2)=S(j1,j2)+SMul2*ic*(FcE+(Fc(i1)+Fc(i2))/3.)
                if(Level.ne.NLevel) then
                  B(i1)=B(i1)-alf*S(j1,j2)*RMassY(i2)
                else
                  If (lOrt) Then
                    call Find(i1,i2,kk,NumNp,MBandD,IAD,IADN)
                    ib=kk
                  Else
                     iB=MBand+i2-i1
                  Endif
                     A(iB,i1)=A(iB,i1)+epsi*S(j1,j2)
                endif
16            continue
17          continue
18          Continue
19        Continue
          Do 20 i=1,NumNP
            If(Level.ne.NLevel) then
              B(i)=B(i)-alf*F(i)
            Else
              if (lOrt) newjjj=IADD(i)
              A(newjjj,i)=A(newjjj,i)+DS(i)/dt
c               A(MBand,i)=A(MBand,i)+DS(i)/dt
              B(i)=B(i)+DS(i)/dt*RMassY(i)-epsi*F(i)
            Endif
20        Continue
21      Continue
C
C  Solve the global matrix equation for transport
C
        If(epsi.lt.0.001) then
          Do 22 i=1,NumNP
          if (lOrt) newjjj = IADD(i)
             B(i)=B(i)/A(newjjj,i)
c            B(i)=B(i)/A(MBand,i)
22        Continue
        Else if(lOrt) then
cccz
c             print*,time,',','Root Diff'
	       call ILU(A,NumNP,MBandD,IAD,IADN,IADD,A1)
             call OrthoMin(A,B1,B,NumNP,MBandD,NumNPD,IAD,
     !	            IADN,IADD,A1,VRV,
     !                RES,RQI,RQ,QQ,QI,RQIDOT,ECNVRG,RCNVRG,ACNVRG,4,
     !                MNorth,MaxItO)
        Else
C
C  Solve the system of linear equations
C
          N1=NumNP-1
          Do 212 k=1,N1
            P=1./A(MBand,k)
            kk=k+1
            kc=MBand
            Do 211 i=kk,NumNP
              kc=kc-1
              If(kc.le.0) goto 212
              C=-P*A(kc,i)
              A(kc,i)=C
              ii=kc+1
              L=kc+MBand-1
              Do 211 j=ii,L
                jj=j+MBand-kc
                A(j,i)=A(j,i)+C*A(jj,k)
211         Continue
212       Continue
          Do  214 i=2,NumNP
            jj=MBand+1-i
            ii=1
            If(jj.le.0) then
              jj=1
              ii=i-MBand+1
            Endif
            Sum=0.
            Do 213 j=jj,MBand-1
              Sum=Sum+A(j,i)*B(ii)
              ii=ii+1
213         Continue
            B(i)=B(i)+Sum
214       Continue
          B(NumNP)=B(NumNP)/A(MBand,NumNP)
          Do 216 k=1,N1
            i=NumNP-k
            jj=i
            m=min0(2*MBand-1,MBand+k)
            Sum=0.
            Do 215 j=MBand+1,m
              jj=jj+1
              Sum=Sum+A(j,i)*B(jj)
215         Continue
            B(i)=(B(i)-Sum)/A(MBand,i)
216       Continue
        Endif
C
        Do i=1,NumNP
          if(lOrt) B(i) = B1(i)
          if (B(i).lt.0) then
                 B(i)=0.0
                 iii=1
            endif
CDT 10/10/11 using the max statement here was causing root system to lose mass. I changed it back.
          RMassY(i)=sngl(B(i))
        Enddo
*

c calculate mature root mass, total root mass and rool lengths
c set up for next step
      Do i=1,NumNPD
        RMean=(RMassY(i)+RMassOld(i))/2.0d0
        RTWT(i)=(RMassM(i)+RMassY(i))
        RMassM(i)=RMassM(i)+ALPY*RMassY(i)*Step
        RMassY(i)=RMassY(i)*(1.0-ALPY*Step)
        RDenM(i)=RMassM(i)/RTWL
        RDenY(i)=RMassY(i)/RTWL
        RDenTotal(i)=RDenM(i)+RDenY(i)

        RMassOld(i)=RMassY(i)

c        VxROld(i)=VxR(i)
c        VzROld(i)=VzR(i)
      Enddo

      TotalRootWeight=0.0D0
      do n=1,NumNP
          TotalRootWeight=TotalRootWeight
     &      +nodeArea(n)*(RMassM(n)+RMassY(n))
      enddo
*
      Return
10    Call errmes(im,il)
      Return
      End
*

      Subroutine RDisper(NumNP,Dispxx,Dispzz,Dispxz,DMolx,DMolz,
     !   hN,Tm,RD)
C      Real xDist, yDist, Dist
CDT 11/2018 we probably don't need this subroutine. All control of the root movement will be
C  done via carbon allocation, we won't be changing the disp coefficiencts based on soil props.
C basically, all this does now is assign values Dispxx and Dispyy
C we can keep it for now if we do want to adjust the values later in the future.
C I have cleaned it up though.
      Dimension
     !  Dispxx(NumNP),Dispzz(NumNP),Dispxz(NumNP),
     !  hN(NumNP),Tm(NumNP),RD(NumNP)

      Do 11 i=1,NumNP
          Adjust=1.0
c          Alp=RD(i)-3.5+3.5*exp(-RD(i)/3.6)
          Alp=1.0D0
          Dispxx(i)=DMolx*Adjust*Alp*min(f1(hN(i)),f2(Tm(i)))
          Dispzz(i)=DMolz*Adjust*Alp*min(f1(hN(i)),f2(Tm(i)))
          Dispxz(i)=0.
11    continue
      Return
      End
*
!      subroutine RWeFact(NumNP,NumEl,NumElD,x,y,KX,WeTab,VxH,VzH,Dispxx,
!     !                  Dispzz,Dispxz)
!      Dimension x(NumNP),y(NumNP),KX(NumElD,4),VxH(NumNP),VzH(NumNP),
!     !          Dispxx(NumNP),Dispzz(NumNP),Dispxz(NumNP),
!     !          WeTab(3,2*NumElD),Beta(3),List(3)
!      Integer e
!
!C AD 10/12/11 commented out the following line, because it causes Vel (geotrophic velocity) to make the denominator very small
!C      TanH(z)=(exp(z)-exp(-z))/(exp(z)+exp(-z))
!      NumSEl=0
!      Do 13 e=1,NumEl
!        NCorn=4
!        If(KX(e,3).eq.KX(e,4)) NCorn=3
!        Do 12 n=1,NCorn-2
!          NumSEl=NumSEl+1
!          M1=KX(e,1)
!          M2=KX(e,n+1)
!          M3=KX(e,n+2)
!          A=y(M2)-y(M1)
!          B=x(M2)-x(M1)
!          Beta(1)=atan2(A,B)
!          A=y(M3)-y(M2)
!          B=x(M3)-x(M2)
!          Beta(2)=atan2(A,B)
!          A=y(M1)-y(M3)
!          B=x(M1)-x(M3)
!          Beta(3)=atan2(A,B)
!          List(1)=M1
!          List(2)=M2
!          List(3)=M3
!          Do 11 j=1,3
!            k=j-1
!            If(k.eq.0) k=3
!            WeTab(k,NumSEl)=0.
!            M1=List(j)
!            jp1=j+1
!            If(j.eq.3) jp1=1
!            M2=List(jp1)
!            Vxx=(VxH(M1)+VxH(M2))/2.
!            Vzz=(VzH(M1)+VzH(M2))/2.
!            If(abs(Vxx).lt.1.e-30.and.abs(Vzz).lt.1.e-30) goto 11
!            BetaV=atan2(Vzz,Vxx)
!            Delta=abs(BetaV-Beta(j))
!            If(Delta.gt.0.314.and.abs(Delta-3.1416).gt.0.314) goto 11
!            ALeng=sqrt((x(M2)-x(M1))**2+(y(M2)-y(M1))**2)
!            CBeta=cos(Beta(j))
!            SBeta=sin(Beta(j))
!            Val=Vxx*CBeta+Vzz*SBeta      ! velocity along the edge
!            VV=sqrt(Vxx*Vxx+Vzz*Vzz)
!            DLL=(Dispxx(M1)+Dispxx(M2))/2.
!            DLT=(Dispxz(M1)+Dispxz(M2))/2.
!            DTT=(Dispzz(M1)+Dispzz(M2))/2.
!            DAL=abs(DLL*CBeta*CBeta+2.0*CBeta*SBeta*DLT+DTT*SBeta*SBeta)
!            Vel=VAL*ALeng
!            Disp=2.0*DAL
!            aa=11.
!            If(Disp.gt.0.) aa=abs(Vel/Disp)
!            If(Disp.lt.1.e-30.or.abs(Vel).lt.0.001*VV.or.abs(aa).gt.10.)
!     !            then
!              If(abs(Vel).lt.0.001*VV) WeTab(k,NumSEl)=0.0
!              If(Vel.gt.0.001*VV) WeTab(k,NumSEl)=1.0
!              If(Vel.lt.-0.001*VV) WeTab(k,NumSEl)=-1.0
!            Else
!              WeTab(k,NumSEl)=1.0/TanH(Vel/Disp)-Disp/Vel
!            Endif
!11        Continue
!12      Continue
!13    Continue
!      Return
!      End

ext: for

max_stars_repo_path: Soil Source/root_diff_new.for

max_stars_repo_name: ARS-CSGCL-DT/MAIZSIM

max_issues_repo_path: Soil Source/root_diff_new.for

max_issues_repo_name: ARS-CSGCL-DT/MAIZSIM

max_forks_repo_path: Soil Source/root_diff_new.for

max_forks_repo_name: ARS-CSGCL-DT/MAIZSIM


text:
C    @(#)senfac.f	20.3 2/13/96
      subroutine senfac

      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha2.inc'
      include 'ipfinc/amtrx.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/beta2.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/gamma.inc'
      include 'ipfinc/slnopt.inc'
      include 'ipfinc/tran.inc'

      common /txstat/ txstat(MAXLTC)
      integer txstat

      integer senltc

      ikec = 1

C     Compute Jacobian elements for LTC

      do jt = 1,ntota
         jndex(1,jt) = ikec
         kstat = senltc (jt,0)
         l = 1

C        Normalize row

         rh = rowh(l)
         rhin = 1.0/rh
         an = rown(l)*rhin
         rj = rowj(l)
         rl = rowl(l)
         rlin = 1.0/(rl-rj*an)
         amtrx(ikec) = rhin
         amtrx(ikec+1) = rj
         amtrx(ikec+2) = rlin
         ikec = ikec+3
         jndex(2,jt) = ikec
         amtrx(ikec) = an
         ikec = ikec + 1
         do l = 2,lp
            mt = kolum(l)
            ah = rowh(l)*rhin
            an = rown(l)*rhin
            amtrx(ikec) = mt
            amtrx(ikec+1) = ah
            amtrx(ikec+2) = (rowj(l) - rj*ah)*rlin
            amtrx(ikec+3) = an
            amtrx(ikec+4) = (rowl(l) - rj*an)*rlin
            ikec = ikec+5
         enddo
      enddo

C     Add an identity row for slack bus elements

      do kta = ntota+1,ntota+nbslck
         jndex(1,kta) = ikec

C        Normalize row

         amtrx(ikec) = 1.0
         amtrx(ikec+1) = 0.0
         amtrx(ikec+2) = 1.0
         ikec = ikec+3
         jndex(2,kta) = ikec
         amtrx(ikec) = 0.0
         ikec = ikec + 1
      enddo

C     Define Jacobian elements for bus constraints

      do kt = nbslck+1,ntot
         kta = kt + ntota
         jndex(1,kta) = ikec
         call senjac (kt,0)
         mel = 0
         mel = korder(mel)
         do while (mel .gt. 0 .and. (kolum(mel) .lt. kta))
            mt = kolum(mel)

C           Eliminate column MT from working row

            ik = jndex(2,mt)
            ikstop = jndex(1,mt+1)
            xn = amtrx(ik)
            ik = ik + 1
            rh = rowh(mel)
            rn = rown(mel)-rh*xn
            rj = rowj(mel)
            rl = rowl(mel)-rj*xn
            amtrx(ikec) = mt
            amtrx(ikec+1) = rh
            amtrx(ikec+2) = rn
            amtrx(ikec+3) = rj
            amtrx(ikec+4) = rl
            ikec = ikec+5
            krw = mel

C           Perform Row MT elimination

            do while (ik .lt. ikstop)
               ml = amtrx(ik)
               if (ml .gt. max) then
                  max = ml
                  ko = lp
                  lp = mend
               else
                  do while (kolum(krw) .lt. ml)
                     ko = krw
                     krw = korder(krw)
                  enddo
               endif
               if (kolum(krw) .eq. ml) then
                  mlc = krw
               else
                  korder(mend) = korder(ko)
                  kolum(mend) = ml
                  korder(ko) = mend
                  mlc = mend
                  ko = mend
                  mend = mend+1
                  rowh(mlc) = 0
                  rown(mlc) = 0
                  rowj(mlc) = 0
                  rowl(mlc) = 0
               endif
               xh = amtrx(ik+1)
               xj = amtrx(ik+2)
               rowh(mlc) = rowh(mlc) - rh*xh - rn*xj
               rowj(mlc) = rowj(mlc) - rj*xh - rl*xj
               xn = amtrx(ik+3)
               xl = amtrx(ik+4)
               rown(mlc) = rown(mlc) - rh*xn - rn*xl
               rowl(mlc) = rowl(mlc) - rj*xn - rl*xl
               ik = ik+5
            enddo
            mel = korder(mel)
         enddo
         if (mel .eq. 0) then

C           Fatal error - no residual diagonal element

            call erexit

         else if (kolum(mel) .gt. kta) then

C           Fatal error - no residual diagonal element

            call erexit

         else

C           Normalize row

            rh = rowh(mel)
            rhin = 1.0/rh
            an = rown(mel)*rhin
            rj = rowj(mel)
            rl = rowl(mel)
            rlin = 1.0/(rl-rj*an)
            amtrx(ikec) = rhin
            amtrx(ikec+1) = rj
            amtrx(ikec+2) = rlin
            ikec = ikec+3
            jndex(2,kta) = ikec
            amtrx(ikec) = an
            ikec = ikec + 1
            mel = korder(mel)
            do while (mel .gt. 0)
               mt = kolum(mel)
               ah = rowh(mel)*rhin
               an = rown(mel)*rhin
               amtrx(ikec) = mt
               amtrx(ikec+1) = ah
               amtrx(ikec+2) = (rowj(mel) - rj*ah)*rlin
               amtrx(ikec+3) = an
               amtrx(ikec+4) = (rowl(mel) - rj*an)*rlin
               ikec = ikec+5
               mel = korder(mel)
            enddo
         endif
      enddo
      jndex(1,ntotx) = ikec

C     End of bus loop

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


subroutine s
  !ERROR: Declaration of 'x' conflicts with its use as internal procedure
  real :: x
contains
  subroutine x
  end
end

module m
  !ERROR: Declaration of 'x' conflicts with its use as module procedure
  real :: x
contains
  subroutine x
  end
end


      PROGRAM LA_GGEV_EXAMPLE

!  -- LAPACK95 EXAMPLE DRIVER ROUTINE (VERSION 1.0) --
!     UNI-C, DENMARK
!     DECEMBER, 1999
!
!  .. "Use Statements" ..
      USE LA_PRECISION, ONLY: WP => SP
      USE F95_LAPACK, ONLY: LA_GGEV
!  .. "Implicit Statement" ..
      IMPLICIT NONE
!  .. "Local Scalars" ..
      INTEGER :: I, J, INFO, N
!  .. "Local Arrays" ..
      COMPLEX(WP), ALLOCATABLE :: A(:,:), AA(:,:), B(:,:), BB(:,:), ALPHA(:)
      COMPLEX(WP), ALLOCATABLE :: VL(:,:), VR(:,:), BETA(:)
      REAL(WP), ALLOCATABLE :: AR(:,:), BR(:,:), ALPHAR(:), ALPHAI(:), BETAR(:), &
                               VLR(:,:), VRR(:,:)
!  .. "Executable Statements" ..
      WRITE (*,*) 'LA_GGEV Example Program Results'
      N = 5
      ALLOCATE( A(N,N), AA(N,N), B(N,N), BB(N,N), ALPHA(N), BETA(N), VL(N,N), VR(N,N) )
      ALLOCATE( AR(N,N), BR(N,N), ALPHAR(N), ALPHAI(N), BETAR(N), VLR(N,N), VRR(N,N) )

      WRITE (*,*) 'Example 1, Real Example'
      OPEN(UNIT=21,FILE='gges.ma',STATUS='UNKNOWN')
      DO J=1,N
         DO I=1,N
            READ(21,*) AR(I,J)
         ENDDO
      ENDDO
      CLOSE(21)
      WRITE(*,*)'Matrix AR : '
      DO I=1,N;
         WRITE(*,"(5(I3,1X))") INT(AR(I,:));
      ENDDO

      OPEN(UNIT=21,FILE='gges.mb',STATUS='UNKNOWN')
      DO J=1,N
         DO I=1,N
            READ(21,*) BR(I,J)
         ENDDO
      ENDDO
      CLOSE(21)
      WRITE(*,*)'Matrix BR : '
      DO I=1,N;
         WRITE(*,"(5(I3,1X))") INT(BR(I,:));
      ENDDO
      WRITE(*,*)
      WRITE(*,*) 'CALL LA_GGEV( A, B, ALPHAR, ALPHAI, BETAR, VLR, VRR )'
      CALL LA_GGEV( AR, BR, ALPHAR, ALPHAI, BETAR, VLR, VRR )

      WRITE(*,*); WRITE(*,*)'ALPHAR on exit : '; WRITE(*,*) ALPHAR(1:N)
      WRITE(*,*); WRITE(*,*)'ALPHAI on exit : '; WRITE(*,*) ALPHAI(1:N)
      WRITE(*,*); WRITE(*,*)'BETA on exit : '; WRITE(*,*) BETAR(1:N)
      WRITE(*,*)'Array VL:'; DO I =1,N; WRITE(*,*)I, VLR(I,1:N); ENDDO
      WRITE(*,*)'Array VR:'; DO I =1,N; WRITE(*,*)I, VRR(I,1:N); ENDDO

      WRITE(*,*)
      WRITE(*,*)' Generalized eigenvalues : '
      DO I=1,N
         WRITE(*,*) '(',ALPHAR(I)/BETAR(I),',',ALPHAI(I)/BETAR(I),')'
      ENDDO

      WRITE (*,*) 'Example 2, Complex Example'
      OPEN(UNIT=21,FILE='ggev.ma',STATUS='UNKNOWN')
      DO J=1,N
         DO I=1,N
            READ(21,*) A(I,J)
         ENDDO
      ENDDO
      CLOSE(21)

      AA=A

      WRITE(*,*)'Matrix A : '
      DO I=1,N
         WRITE(*,"(5('('(I3,1X,',',I3)')',1X,1X))") INT(A(I,:)), INT(AIMAG(A(I,:)))
      ENDDO

      OPEN(UNIT=21,FILE='ggev.mb',STATUS='UNKNOWN')
      DO J=1,N
         DO I=1,N
            READ(21,*) B(I,J)
         ENDDO
      ENDDO
      CLOSE(21)

      BB=B

      WRITE(*,*)'Matrix B : '
      DO I=1,N
         WRITE(*,"(5('('(I3,1X,',',I3)')',1X,1X))") INT(B(I,:)), INT(AIMAG(B(I,:)))
      ENDDO
      WRITE(*,*)
      WRITE(*,*) 'CALL LA_GGEV( A, B, ALPHA, BETA, INFO=INFO )'
      CALL LA_GGEV( A, B, ALPHA, BETA,  INFO=INFO )

      WRITE(*,*)
      WRITE(*,*)'ALPHA on exit : '
      DO I=1,N
         WRITE(*,*) ALPHA(I)
      ENDDO

      WRITE(*,*)
      WRITE(*,*)'BETA on exit : '
      DO I=1,N
         WRITE(*,*) BETA(I)
      ENDDO
      WRITE(*,*)
      WRITE(*,*)'INFO = ', INFO

      WRITE(*,*)
      WRITE(*,*)' Generalized eigenvalues : '
      DO I=1,N
         WRITE(*,*)  ALPHA(I)/BETA(I)
      ENDDO

      WRITE(*,*)
!     WRITE(*,*)' * EXAMPLE 2 * '
      WRITE(*,*) 'CALL LA_GGEV( A, B, ALPHAR, ALPHAI, BETA, VL, VR )'
      CALL LA_GGEV( A, B, ALPHA, BETA, VL, VR )

      WRITE(*,*)
      WRITE(*,*)'Matrix VL on exit : '
      DO I=1,N;
         WRITE(*,"(5('('(F8.5,1X,',',F8.5)')'))") REAL(VL(I,:)), AIMAG(VL(I,:))
      ENDDO

      WRITE(*,*)
      WRITE(*,*)'Matrix VR on exit : '
      DO I=1,N;
         WRITE(*,"(5('('(F8.5,1X,',',F8.5)')'))") REAL(VR(I,:)), AIMAG(VR(I,:))
      ENDDO

      END PROGRAM LA_GGEV_EXAMPLE


text:
while_datatypes/main.cpp
while_datatypes/stimulus.cpp
while_datatypes/display.cpp
while_datatypes/while_datatypes.cpp


text:
The HäagenDazs Honey Bee Haven is a UC Davis garden designed to provide nourishment to honeybees.  It was created as the result of a $125,000 gift to the UC Davis Department of Entomology from the HäagenDazs ice cream company. A major goal of the garden is to raise public awareness about Colony Collapse Disorder (CCD), a mysterious phenomenon which may pose a serious threat to bee populations.

The gardens grand opening took place on September 11, 2010. Its currently open to the public all yeararound, whenever the sun is up. Its free to get in.

The garden is located just east of the Harry H. Laidlaw Jr. Honey Bee Research Facility, and its area encompasses the site upon which the Baxter House once stood.

 Photos from the HäagenDazs Honey Bee Haven

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C
C +-----------------------------------------------------------------+
C |                                                                 |
C |                Copyright (C) 1986 by UCAR                       |
C |        University Corporation for Atmospheric Research          |
C |                    All Rights Reserved                          |
C |                                                                 |
C |                 NCARGRAPHICS  Version 1.00                      |
C |                                                                 |
C +-----------------------------------------------------------------+
C
C
C ---------------------------------------------------------------------
C
      SUBROUTINE AGEXAX (IAXS,SVAL,UMIN,UMAX,NICE,QLUA,FUNS,QBTP,BASD,
     +                   BASE,QMJD,QMND,QMNT,QLTD,QLTP,QLED,QLEX,QLFD,
     +                   QLFL,QMIN,QMAX)
C
      DIMENSION SVAL(2)
C
C The routine AGEXAX is used by AGSTUP to examine the parameters which
C determine how a given axis is tick-marked and labelled and to provide
C default values for missing ones.  Its arguments are as follows:
C
C -- IAXS is the number of the axis being drawn - 1, 2, 3, or 4.
C
C -- SVAL is the array of special values.
C
C -- UMIN and UMAX are the minimum and maximum values along the axis, in
C    the user coordinate system.  Rounded values of UMIN and UMAX are
C    returned in QMIN and QMAX if the following argument (NICE) is zero.
C
C -- NICE is a flag indicating whether rounded values of UMIN and UMAX
C    are to be returned (NICE.EQ.0) or not (NICE.NE.0).
C
C -- LLUA and FUNS specify the user-system-to-label-system mapping along
C    the axis.  See the routine AGAXIS for a discussion of them.
C
C -- NBTP, BASD, BASE, and NMJD are used to determine the positioning of
C    major tick marks in the label coordinate system.  NBTP and BASE are
C    described in the routine AGNUMB.  BASD is the desired value of BASE
C    supplied by the user.  If BASD has a null value, BASE is computed
C    by AGEXAX.  NMJD is a user-supplied-or-defaulted parameter giving
C    the approximate number of major ticks (and therefore the number of
C    numeric labels) to be placed on the axis.
C
C -- NMND and NMNT are the desired and actual (to be determined) number
C    of minor ticks per major division.  See discussion in AGAXIS.
C
C -- NLTD, NLTP, NLED, NLEX, NLFD, and NLFL are desired and actual (to
C    be determined) values of the parameters describing the form to be
C    used for numeric labels.  See discussion in AGNUMB.
C
C -- QMIN and QMAX are rounded values of UMIN and UMAX, returned only if
C    NICE.EQ.0.
C
C The following common block contains AUTOGRAPH variables which are
C not control parameters.  The only one used here is SMRL, which is a
C (machine-dependent) small real which, when added to a number in the
C range (1,10), will round it upward without seriously affecting the
C leading significant digits.  The object of this is to get rid of
C strings of nines.
C
      COMMON /AGORIP/ SMRL , ISLD , MWCL,MWCM,MWCE,MDLA,MWCD,MWDQ ,
     +                INIF
C
C The arrays BASP and NMNP specify possible default values for BASE and
C NMNT when NBTP.EQ.1.
C
      DIMENSION BASP(5),NMNP(5)
C
      DATA BASP(1) / 10. / , NMNP(1) / 1 / ,
     *     BASP(2) /  5. / , NMNP(2) / 4 / ,
     *     BASP(3) /  2. / , NMNP(3) / 1 / ,
     *     BASP(4) /  1. / , NMNP(4) / 1 / ,
     *     BASP(5) /  .5 / , NMNP(5) / 4 /
C
C If the parameter NBTP is zero, tick marks and labels are suppressed.
C
      NBTP=IFIX(QBTP)
      IF (NBTP.EQ.0) RETURN
C
C Unpack integer values from floating-point arguments.
C
      LLUA=IFIX(QLUA)
      NMJD=IFIX(QMJD)
      IF (QMND.NE.SVAL(1).AND.QMND.NE.SVAL(2)) NMND=IFIX(QMND)
      NMNT=0
      IF (QLTD.NE.SVAL(1).AND.QLTD.NE.SVAL(2)) NLTD=IFIX(QLTD)
      NLTP=0
      IF (QLED.NE.SVAL(1).AND.QLED.NE.SVAL(2)) NLED=IFIX(QLED)
      NLEX=0
      IF (QLFD.NE.SVAL(1).AND.QLFD.NE.SVAL(2)) NLFD=IFIX(QLFD)
      NLFL=0
C
C Compute label-coordinate-system values at the ends of the axis.
C
      CALL AGUTOL (IAXS,FUNS,1,UMIN,VMIN)
      CALL AGUTOL (IAXS,FUNS,1,UMAX,VMAX)
C
C Error if the label-coordinate-system values are equal.
C
      IF (VMIN.EQ.VMAX) GO TO 901
C
C If a special value is specified for the parameter BASD, AGEXAX must
C pick a value for the parameter BASE.
C
      IF (BASD.EQ.SVAL(1).OR.BASD.EQ.SVAL(2)) GO TO 101
C
C The user has specified a value for the parameter BASE.  If that value
C is less than or equal to zero, tick marks and labels are suppressed.
C
      BASE=AMAX1(0.,BASD)
      IF (BASE.EQ.0.) RETURN
      NMNT=0
      GO TO 108
C
C Pick a value for the parameter BASE, depending on the number type.
C
  101 GO TO (102,105,106) , NBTP
C
C Major ticks and labels are at numbers of the form (-) BASE * EXMU.
C
  102 NMJD=MAX0(0,NMJD)
C
C Compute an approximate value for BASE.
C
      FTMP=ABS(VMAX-VMIN)/FLOAT(NMJD+1)
C
C Reduce the approximate value to the form FTMP * 10 ** ITMP.
C
      ASSIGN 103 TO JMP1
      GO TO 200
C
C Pick a reasonable value for BASE (1., 2., OR 5. * 10**ITMP).
C
  103 DO 104 I=1,5
        IF (FTMP.LT.BASP(I)) GO TO 104
        BASE=BASP(I)*SNGL(10.D0**ITMP)
        NMNT=NMNP(I)
        GO TO 107
  104 CONTINUE
C
C Major ticks and labels are at numbers of the form (-) BASE * 10**EXMU.
C
  105 BASE=1.
      NMNT=8
      GO TO 107
C
C Major ticks and labels are at numbers of the form (-) BASE**EXMU.
C
  106 BASE=10.
      NMNT=8
C
  107 IF (BASD.EQ.SVAL(2)) BASD=BASE
C
  108 IF (QMND.NE.SVAL(1).AND.QMND.NE.SVAL(2)) NMNT=MAX0(0,NMND)
      IF (QMND.EQ.SVAL(2)) QMND=FLOAT(NMNT)
C
C If the user wants nice values at the axis ends, reset UMIN and UMAX.
C
      IF (NICE.NE.0) GO TO 115
C
      LOOP=0
C
      WMIN=VMIN
      WMAX=VMAX
C
      GO TO (109,110,112) , NBTP
C
  109 EMIN=VMIN/BASE+.5+SIGN(.5,VMIN-VMAX)
      EMIN=EMIN-.5+SIGN(.5,EMIN)-SIGN(SMRL*EMIN,VMIN-VMAX)
      WMIN=BASE*(EMIN-AMOD(EMIN,1.))
      EMAX=VMAX/BASE+.5+SIGN(.5,VMAX-VMIN)
      EMAX=EMAX-.5+SIGN(.5,EMAX)-SIGN(SMRL*EMAX,VMAX-VMIN)
      WMAX=BASE*(EMAX-AMOD(EMAX,1.))
      GO TO 114
C
  110 IF (VMIN.EQ.0.) GO TO 111
      EMIN=ALOG10(ABS(VMIN)/BASE)+.5+SIGN(.5,VMIN*(VMIN-VMAX))
      EMIN=EMIN-.5+SIGN(.5,EMIN)-SIGN(SMRL*EMIN,VMIN*(VMIN-VMAX))
      WMIN=SIGN(BASE,VMIN)*10.**(EMIN-AMOD(EMIN,1.))
  111 IF (VMAX.EQ.0.) GO TO 114
      EMAX=ALOG10(ABS(VMAX)/BASE)+.5+SIGN(.5,VMAX*(VMAX-VMIN))
      EMAX=EMAX-.5+SIGN(.5,EMAX)-SIGN(SMRL*EMAX,VMAX*(VMAX-VMIN))
      WMAX=SIGN(BASE,VMAX)*10.**(EMAX-AMOD(EMAX,1.))
      GO TO 114
C
  112 IF (BASE.EQ.1.) GO TO 115
      IF (VMIN.EQ.0.) GO TO 113
      EMIN=ALOG10(ABS(VMIN))/ALOG10(BASE)+.5+SIGN(.5,VMIN*(VMIN-VMAX))
      EMIN=EMIN-.5+SIGN(.5,EMIN)-SIGN(SMRL*EMIN,VMIN*(VMIN-VMAX))
      WMIN=SIGN(1.,VMIN)*BASE**(EMIN-AMOD(EMIN,1.))
  113 IF (VMAX.EQ.0.) GO TO 114
      EMAX=ALOG10(ABS(VMAX))/ALOG10(BASE)+.5+SIGN(.5,VMAX*(VMAX-VMIN))
      EMAX=EMAX-.5+SIGN(.5,EMAX)-SIGN(SMRL*EMAX,VMAX*(VMAX-VMIN))
      WMAX=SIGN(1.,VMAX)*BASE**(EMAX-AMOD(EMAX,1.))
C
C Re-compute the user-coordinate-system minimum and maximum values.
C
  114 CALL AGUTOL (IAXS,FUNS,-1,WMIN,QMIN)
      CALL AGUTOL (IAXS,FUNS,-1,WMAX,QMAX)
C
C Test for problems with nice values chosen.
C
      IF (QMIN.LT.QMAX) GO TO 140
      IF (QMIN.GT.QMAX) GO TO 901
C
C We have a pathological case - user values are clustered very close to
C a label position.  See what can be done about it.
C
      LOOP=LOOP+1
      IF (LOOP.GT.1) GO TO 901
C
      GO TO (137,138,139) , NBTP
C
  137 VMIN=VMIN+SIGN(BASE,VMIN-VMAX)
      VMAX=VMAX+SIGN(BASE,VMAX-VMIN)
      GO TO 109
C
  138 VMIN=VMIN*10.**SIGN(1.,VMIN*(VMIN-VMAX))
      VMAX=VMAX*10.**SIGN(1.,VMAX*(VMAX-VMIN))
      GO TO 110
C
  139 VMIN=VMIN*BASE**SIGN(1.,VMIN*(VMIN-VMAX))
      VMAX=VMAX*BASE**SIGN(1.,VMAX*(VMAX-VMIN))
      GO TO 112
C
  140 VMIN=WMIN
      VMAX=WMAX
C
C Now we examine the parameters defining the appearance of the numeric
C labels.  If the numeric-label type is zero, there is no more to do.
C
  115 IF (QLTD.EQ.SVAL(1).OR.QLTD.EQ.SVAL(2)) GO TO 116
      NLTP=MAX0(0,MIN0(3,NLTD))
      IF (NLTP.EQ.0) GO TO 136
C
C The numeric-label type (NLTP) is specified.  If both the numeric-label
C exponent and numeric-label fraction-length are also specified, quit.
C
      NLEX=NLED
      NLFL=NLFD
      IF (QLED.NE.SVAL(1).AND.QLED.NE.SVAL(2).AND.
     +    QLFD.NE.SVAL(1).AND.QLFD.NE.SVAL(2)       ) GO TO 136
      GO TO 117
C
C We must pick a value for the numeric-label type.  Start with the dummy
C value 4 so as to jump to the proper piece of code.
C
  116 NLTP=4
C
C Reduce the value of BASE to the form RBSE * 10**KBSE, where RBSE is
C in the range (1,10) and KBSE is an integer.
C
  117 FTMP=BASE
      ASSIGN 118 TO JMP1
      GO TO 200
C
  118 RBSE=FTMP
      KBSE=ITMP
C
C Compute LBSE = the number of significant digits in RBSE.
C
      ASSIGN 119 TO JMP2
      GO TO 300
C
  119 LBSE=1+ITMP
C
C Jump depending on the value of the numeric-label type.
C
      GO TO (120,128,131,132) , NLTP
C
C Scientific notation is to be used.  Estimate the number of significant
C digits that are likely to be required, depending on the number type.
C
  120 GO TO (121,123,124) , NBTP
C
  121 FTMP=AMAX1(ABS(VMIN),ABS(VMAX))/BASE
      ASSIGN 122 TO JMP1
      GO TO 200
C
  122 NSIG=MAX0(1,ITMP+1+LBSE)
      GO TO 125
C
  123 NSIG=LBSE
      GO TO 125
C
  124 NSIG=10
C
C NLEX + NLFL should be equal to NSIG.  Make that the case.
C
  125 IF (QLED.NE.SVAL(1).AND.QLED.NE.SVAL(2)) GO TO 127
      IF (QLFD.EQ.SVAL(1).OR. QLFD.EQ.SVAL(2)) GO TO 126
      NLEX=NSIG-MAX0(0,NLFL)
      GO TO 135
  126 NLEX=1
  127 NLFL=NSIG-NLEX
      IF (NLFL.LE.0) NLFL=-1
      GO TO 135
C
C Exponential notation is to be used.  Compute the exponent NEXP such
C that BASE / 10**NEXP is an integer.
C
  128 NEXP=KBSE-LBSE+1
C
C NLEX - NLFL should be equal to NEXP.  Make that the case.  (Note that,
C if NBTP is 3, NLEX is forced to zero.)
C
      IF (NBTP.EQ.3) NLEX=0
C
      IF (QLFD.NE.SVAL(1).AND.QLFD.NE.SVAL(2)) GO TO 129
      IF (QLED.NE.SVAL(1).AND.QLED.NE.SVAL(2)) GO TO 130
      NLFL=-1
  129 NLEX=MAX0(0,NLFL)+NEXP
      GO TO 135
  130 NLFL=NLEX-NEXP
      IF (NLFL.LE.0) NLFL=-1
      GO TO 135
C
C No-exponent notation is to be used.  NLFL is the only parameter we
C need to worry about.  If it is already set, quit.
C
  131 IF (QLFD.NE.SVAL(1).AND.QLFD.NE.SVAL(2)) GO TO 136
C
C Set NLFL to the actual number of digits in the fractional portion of
C BASE.
C
      NLFL=LBSE-KBSE-1
      IF (NLFL.LE.0) NLFL=-1
      GO TO 135
C
C We must pick a value for the numeric-label type, depending on the
C number type.
C
  132 GO TO (133,134,134) , NBTP
C
C Nunbers are of the form (-) BASE * EXMU.  Use labels with no exponent
C unless the use of an exponent would result in shorter labels.
C
  133 IF (MAX0(KBSE+1-LBSE,-KBSE-1).GT.4) GO TO 134
      NLTP=3
      NLFL=LBSE-KBSE-1
      IF (NLFL.LE.0) NLFL=-1
      GO TO 135
C
C Exponential notation is used.
C
  134 NLTP=2
      NLEX=KBSE-LBSE+1
      NLFL=-1
C
C Back-store the computed parameters, if requested, and return.
C
  135 IF (QLTD.EQ.SVAL(2)) QLTD=FLOAT(NLTP)
      IF (QLED.EQ.SVAL(2)) QLED=FLOAT(NLEX)
      IF (QLFD.EQ.SVAL(2)) QLFD=FLOAT(NLFL)
C
C Pack up integer values to floating-point arguments and return.
C
  136 QMNT=FLOAT(NMNT)
      QLTP=FLOAT(NLTP)
      QLEX=FLOAT(NLEX)
      QLFL=FLOAT(NLFL)
      RETURN
C
C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
C
C This internal procedure reduces the number (FTMP) to the range (1,10),
C returning (FTMP) and (ITMP) such that (FTMP) * 10**(ITMP) is equal to
C the original value of (FTMP).  (FTMP) must be positive.
C
  200 FTM1=ALOG10(FTMP+SMRL*FTMP)
      IF (FTM1.LT.0.) FTM1=FTM1-1.
      ITMP=IFIX(FTM1)
      FTMP=AMAX1(1.,FTMP*SNGL(10.D0**(-ITMP)))
      GO TO JMP1 , (103,118,122)
C
C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
C
C This internal procedure counts the number of digits in the fractional
C portion of (FTMP), returning the count as the value of (ITMP).
C
  300 FTM1=AMOD(FTMP+SMRL*FTMP,1.)
      FTM2=10.*SMRL*FTMP
      ITMP=0
C
  301 IF (FTM1.LT.FTM2) GO TO 302
      ITMP=ITMP+1
      IF (ITMP.GE.10) GO TO 302
      FTM1=AMOD(10.*FTM1,1.)
      FTM2=10.*FTM2
      GO TO 301
C
  302 GO TO JMP2 , (119)
C
C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
C
C Error exit.
C
C +NOAO - Comment out FTN write and format statement, SETER is okay.
C
  901 CONTINUE
C 901 WRITE (I1MACH(4),9001) IAXS
      CALL SETER ('AGEXAX (CALLED BY AGSTUP) - USER-SYSTEM-TO-LABEL-SYST
     +EM MAPPING IS NOT MONOTONIC',1,2)
C
C Formats.
C
C9001 FORMAT ('0PROBLEM WITH AXIS NUMBER',I2,
C    +        ' (1, 2, 3, AND 4 IMPLY LEFT, RIGHT, BOTTOM, AND TOP)')
C
C -NOAO
      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
program files
use rm
implicit none

character(len=100)::file="write.dat"

integer::nrows=10,ncols=10,iniseed=3,average=1
real::npos=0.5,nneg=0.5
integer,allocatable::matriz(:,:)
!Index to loop the matrix
integer::i,j,k,IDfile1 = 35
real::rvec(1)

allocate(matriz(ncols,nrows))
matriz(:,:) = 0
do j=1,nrows
	do i = 1,ncols
		call ranmar(rvec,1,iniseed*average)
		!Initial seed set to iniseed*average so that the programs
		!runs with a different initial configuration depending
		!on the average value
		if (rvec(1) < npos) then
			matriz(i,j) = 1
		else
			matriz (i,j) = -1
		end if
	end do
end do

open (UNIT=IDfile1, FILE=file,ACTION='write',STATUS='replace')
!close(IDfile1)
!open (UNIT=IDfile1, FILE=file,ACTION='write',STATUS='append')
do k = 1,2
	do j = 1,nrows
		do i = 1,ncols

			if(j== nrows .AND. i== ncols) then
				write(35,'(I2.1,/)',advance='no') matriz(i,j)
				print *, i,j
			!else if(j == 1 .AND. i==1) then
			!	write(35,'(I2.1)') matriz(i,j)
			!	write(35,'(A)',advance='no') ","
			!	print *, i,j
			else
				write(35,'(I2.1)',advance='no') matriz(i,j)
				write(35,'(A)',advance='no') ","
			end if
		end do
	end do
end do
close(IDfile1)
end program files

max_stars_repo_name: mattborghi/Thesis

max_issues_repo_name: mattborghi/Thesis

max_forks_repo_name: mattborghi/Thesis


text:
!-----------------------------------------------------------------------------!
! \file   kba/Block_Tracer.f90
! \author Thomas M. Evans
! \date   Thur Jan 24 14:00 2008
! \brief  FORTRAN kernel for tracing a ray through a block of mesh.
! \note   Copyright (C) 2008 Oak Ridge National Laboratory, UT-Battelle, LLC
!-----------------------------------------------------------------------------!
! $Id: Ray_Tracer_Kernel.f90,v 1.1 2008/02/07 21:27:48 9te Exp $
!-----------------------------------------------------------------------------!

module Ray_Tracer

  use Parameters, only : DOUBLE
  implicit none
  public

  ! Distance-to-surface arrays used by the ray tracer.
  real(DOUBLE), dimension(:), allocatable, save :: dx, dy, dz

  ! Ray direction.
  real(DOUBLE), dimension(3), save :: omega

contains

  !---------------------------------------------------------------------------!
  ! Calculate distance to boundary arrays
  !---------------------------------------------------------------------------!

  subroutine calc_distances(start, edges, r, p, angle, d)

    implicit none

    ! >>> INPUT/OUTPUT DATA

    ! Starting cell position in edge array.
    integer, intent(in) :: start

    ! Edges in the chosen direction.
    real(DOUBLE), dimension(:), intent(in) :: edges

    ! Starting point position in the chosen direction.
    real(DOUBLE), intent(in) :: r

    ! Target (ending) point position in the chosen direction.
    real(DOUBLE), intent(in) :: p

    ! Inverse of the angle-of-flight.
    real(DOUBLE), intent(in) :: angle

    ! Edges in the chosen direction.
    real(DOUBLE), dimension(:), intent(out) :: d

    ! >>> LOCAL DATA

    ! Indices.
    integer :: n, index

    ! Number of edges.
    integer :: length

    ! Inverse angle.
    real(DOUBLE) :: inv_angle

    ! >>> BODY

    ! intialize d
    d = 1.0e14

    ! calculate the size of the edges array
    length = size(edges)

    ! initialize starting index
    n     = start
    index = 1

    ! step through edges in the positive direction
    if (angle > 0.0) then

       ! calculate inverse of angle
       inv_angle = 1.0_DOUBLE / angle

       ! continue until we hit the point or we get to the edge of the mesh
       do while (edges(n) < p .and. n < length)

          ! assign distance to boundary
          d(index) = (edges(n + 1) - r) * inv_angle

          ! update the indices
          n     = n + 1
          index = index + 1

       end do

    else if (angle < 0.0) then

       ! calculate inverse of angle
       inv_angle = 1.0_DOUBLE / angle

       ! continue until we hit the point or we get to the edge of the mesh
       do while (edges(n + 1) > p .and. n > 0)

          ! assign distance to boundary
          d(index) = (edges(n) - r) * inv_angle

          ! update the indices
          n     = n - 1
          index = index + 1

       end do

    end if

  end subroutine calc_distances

  !---------------------------------------------------------------------------!
  ! CONSTRUCT
  !---------------------------------------------------------------------------!

  subroutine construct_ray_tracer

    use Parameters, only : im, jm, km
    implicit none

    if (.not. allocated(dx)) allocate(dx(im + 1))
    if (.not. allocated(dy)) allocate(dy(jm + 1))
    if (.not. allocated(dz)) allocate(dz(km + 1))

  end subroutine construct_ray_tracer

  !---------------------------------------------------------------------------!
  ! DESTRUCT
  !---------------------------------------------------------------------------!

  subroutine destruct_ray_tracer

    implicit none

    if (allocated(dx)) deallocate(dx)
    if (allocated(dy)) deallocate(dy)
    if (allocated(dz)) deallocate(dz)

  end subroutine destruct_ray_tracer

end module Ray_Tracer

!-----------------------------------------------------------------------------!
! Setup routines for ray-tracing
!-----------------------------------------------------------------------------!

subroutine setup_ray_tracing

  use Ray_Tracer, only : construct_ray_tracer

  call construct_ray_tracer

end subroutine setup_ray_tracing

!-----------------------------------------------------------------------------!

subroutine cleanup_ray_tracing

  use Ray_Tracer, only : destruct_ray_tracer

  call destruct_ray_tracer

end subroutine cleanup_ray_tracing

!-----------------------------------------------------------------------------!
! Ray-tracing through a block mesh
!-----------------------------------------------------------------------------!

! Trace through a mesh block.

subroutine block_tracer(x, y, z, sigma, r, p, ijk, tau, terminator)

  use Parameters, only : SHORT, DOUBLE, im, jm, km, I, J, K, one
  use Ray_Tracer, only : dx, dy, dz, omega, calc_distances
  implicit none

  ! >>> INPUT DATA

  ! Edge arrays.
  real(DOUBLE), dimension(im + 1), intent(in) :: x
  real(DOUBLE), dimension(jm + 1), intent(in) :: y
  real(DOUBLE), dimension(km + 1), intent(in) :: z

  ! Total cross sections in each cell for all groups in (1/cm).
  real(DOUBLE), dimension(im * jm * km), intent(in) :: sigma

  ! Ray target (ending) point.
  real(DOUBLE), dimension(3), intent(in) :: p

  ! >>> IN/OUT DATA

  ! (i, j, k) position of ray in mesh.
  integer, dimension(3), intent(inout) :: ijk

  ! Ray position.
  real(DOUBLE), dimension(3), intent(inout) :: r

  ! Tau = sum_{i} sigma_i * l_i where l is the length of each step.
  real(DOUBLE), intent(inout) :: tau

  ! >>> OUT DATA

  ! Termination condition.
  !   0 - hit target point
  !   1 - exited low-X face of block
  !   2 - exited high-X face of block
  !   3 - exited low-Y face of block
  !   4 - exited high-Y face of block
  !   5 - exited low-Z face of block
  !   6 - exited high-Z face of block
  !   7 - negative distance-to-boundary
  integer, intent(out) :: terminator

  ! >>> LOCAL DATA

  ! Distance to the source point: norm = 1 / ||(r - rp)||.
  real(DOUBLE) :: dtarget, dtravel, dstep, d

  ! Boolean flag for ray-tracing.
  logical :: done

  ! Indices into distance arrays.
  integer :: ii, jj, kk

  ! Canonical cell index.
  integer :: cell, step_i, step_j, step_k

  ! >>> BODY

  ! add 1 to ijk to adjust to FORTRAN-style indexing
  ijk  = ijk + 1

  ! calculate the distance to the target and the direction of the ray
  dtarget  = sqrt(dot_product(r - p, r - p))
  omega    = (p - r) / dtarget

  ! calculate the distances to each surface
  call calc_distances(ijk(I), x, r(I), p(I), omega(I), dx)
  call calc_distances(ijk(J), y, r(J), p(J), omega(J), dy)
  call calc_distances(ijk(K), z, r(K), p(K), omega(K), dz)

  ! initialize indices into distance arrays
  ii = 1
  jj = 1
  kk = 1

  ! steps
  step_i = int(sign(one, omega(I)))
  step_j = int(sign(one, omega(J)))
  step_k = int(sign(one, omega(K)))

  ! initialize loop
  done    = .false.
  dtravel = 0.0_DOUBLE

  ! ray-trace to point or until we leave the block
  do while (done .eqv. .false.)

     ! calculate the canonical cell index
     cell = ijk(I) + (ijk(J) - 1) * im + (ijk(K) - 1) * im * jm

     ! get the distance
     d = min(dx(ii), min(dy(jj), dz(kk)))

     ! calculate the step
     dstep = d - dtravel

     ! push to next boundary or to target point
     if (dstep < dtarget) then

        ! subtract the distance from the target and add to the travel
        dtarget = dtarget - dstep
        dtravel = d

        ! add up tau
        tau = tau + sigma(cell) * dstep

        ! step to next face
        if (d == dx(ii)) then

           ii     = ii + 1
           ijk(I) = ijk(I) + step_i

           ! terminating conditions
           if (ijk(I) == 0) then
              terminator = 1
              done       = .true.
           else if (ijk(I) == im + 1) then
              terminator = 2
              done       = .true.
           end if

        else if (d == dy(jj)) then

           jj     = jj + 1
           ijk(J) = ijk(J) + step_j

           ! terminating conditions
           if (ijk(J) == 0) then
              terminator = 3
              done       = .true.
           else if (ijk(J) == jm + 1) then
              terminator = 4
              done       = .true.
           end if

        else if (d == dz(kk)) then

           kk     = kk + 1
           ijk(K) = ijk(K) + step_k

           ! terminating conditions
           if (ijk(K) == 0) then
              terminator = 5
              done       = .true.
           else if (ijk(K) == km + 1) then
              terminator = 6
              done       = .true.
           end if

        end if

     else

        ! add to the travel
        dtravel = dtravel + dtarget

        ! add up tau
        tau = tau + sigma(cell) * dtarget

        ! we are finished
        done       = .true.
        terminator = 0

     end if

  end do

  ! calculate final location of particle
  r = r + dtravel * omega

  ! subtract 1 to ijk to adjust to C-style indexing
  ijk = ijk - 1

end subroutine block_tracer

!-----------------------------------------------------------------------------!
!                          end of Block_Tracer.f90
!-----------------------------------------------------------------------------!

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


           program demo_get_tmp
           use M_io, only : get_tmp, uniq
           implicit none
           character(len=:),allocatable :: answer
              answer=get_tmp()
              write(*,*)'result is ',answer
              answer=get_tmp()//uniq('_scratch',create=.false.)
              write(*,*)'the file ',answer,' was a good scratch file name, at least a moment ago'
           end program demo_get_tmp