Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran77 woes- single parameter least squares minimisation stuck!

  1. Jan 26, 2014 #1
    Hey guys! Having issues with a code I'm writing. It'll (eventually) be a four-parameter least squares minimisation, but for the moment I'm just trying to get the single parameter version to work. I've got the first module done but, when I run it it gets caught up in a loop and doesn't actually minimise anything! Good start huh?

    Here's my code, bless me with your insight!

    Code (Text):
          PROGRAM assignment
         
          ! A program designed to fit experiemental data, using the method
          ! of least squares to minimise the associated chi-squared and
          ! obtain the four control parameters A,B,h1 and h2.
          !*****************************************************************

          IMPLICIT NONE
          INTEGER i,t0
          DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
          DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
          DOUBLE PRECISION deltah,Cs(17)

            OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
               DO i=1,17
                  READ(21,*)t(i),Ct(i),eCt(i)
               END DO
            CLOSE(21)

          !Read in data.txt as three one dimensional arrays.
          !*****************************************************************
          !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
          !DO i=1,17
          !   WRITE(21,*)t(i),Ct(i),eCt(i)
          !END DO
          !CLOSE(21)
          !
          !Just to check input file is being read correctly.
          !*****************************************************************
          !****************************Parameters***************************
         
          Aa= 0
          Ba= 0
          h1a= 0
          h2a= 0

          !*****************************************************************
          !**********************Minimising Lamda1 (h1)*********************
          deltah= 0
               DO 10
          chisqa= 0
                  DO 20 i= 1,17
                     h1a= h1a + deltah
                       Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
                         Ctdiff= Ct - Cs
                           c= Ctdiff**2/eCt**2
                             chisqa= chisqa + c(i)
                     deltah= 0.001
     20           END DO
          !Use initial h1 value of 0 to calculate start-point chisq.

          h1b= 0.001
          chisqb= 0
                  DO 30 i= 1,17
                       h1b= h1b + deltah
                         Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
                           Ctdiff(i)= Ct(i) - Cs(i)
                             c= Ctdiff(i)**2/eCt(i)**2
                               chisqb= chisqb + c(i)
     30           END DO
          !First-step h1 used to find competing chisq for comparison.
                     WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
                     WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
          !Prints the two calculated chisq values to screen.
                     dchisq= chisqa - chisqb
                  IF (dchisq.GT.0) THEN
                     h1a= h1b
                     ELSE IF (dchisq.LE.0) THEN
                     deltah= (-deltah*2)/100

                  END IF
               IF (chisqb.LE.6000) EXIT
     10        END DO
               WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
          !*****************************************************************
               
          END PROGRAM assignment

    Thanks for the help.
    Nathan

    EDIT: I should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!


    EDIT: I've read this is the wrong place for this. Sorry, have opened a thread in the correct place, please remove this one.
     
    Last edited: Jan 26, 2014
  2. jcsd
  3. Jan 26, 2014 #2
    Code (Text):

         1  PROGRAM assignment
         2 
         3  ! A program designed to fit experiemental data, using the method
         4  ! of least squares to minimise the associated chi-squared and
         5  ! obtain the four control parameters A,B,h1 and h2.
         6  !*****************************************************************
         7 
         8  IMPLICIT NONE
         9  INTEGER i,t0
        10  DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
        11  DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
        12  DOUBLE PRECISION deltah,Cs(17)
        13 
        14  OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
        15  DO i=1,17
        16  READ(21,*)t(i),Ct(i),eCt(i)
        17  END DO
        18  CLOSE(21)
        19 
        20  !Read in data.txt as three one dimensional arrays.
        21  !*****************************************************************
        22  !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
        23  !DO i=1,17
        24  ! WRITE(21,*)t(i),Ct(i),eCt(i)
        25  !END DO
        26  !CLOSE(21)
        27  !
        28  !Just to check input file is being read correctly.
        29  !*****************************************************************
        30  !****************************Parameters***************************
        31 
        32  Aa= 0
        33  Ba= 0
        34  h1a= 0
        35  h2a= 0
        36 
        37  !*****************************************************************
        38  !**********************Minimising Lamda1 (h1)*********************
        39  deltah= 0
        40  DO 10
        41      chisqa= 0
        42      DO 20 i= 1,17
        43          h1a= h1a + deltah
        44          Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
        45          Ctdiff= Ct - Cs
        46          c= Ctdiff**2/eCt**2
        47          chisqa= chisqa + c(i)
        48          deltah= 0.001
        49      20 END DO
        50     
        51      !Use initial h1 value of 0 to calculate start-point chisq.
        52      h1b= 0.001
        53      chisqb= 0
        54      DO 30 i= 1,17
        55          h1b= h1b + deltah
        56          Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
        57          Ctdiff(i)= Ct(i) - Cs(i)
        58          c= Ctdiff(i)**2/eCt(i)**2
        59          chisqb= chisqb + c(i)
        60      30 END DO
        61 
        62      !First-step h1 used to find competing chisq for comparison.
        63      WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
        64      WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
        65      !Prints the two calculated chisq values to screen.
        66      dchisq= chisqa - chisqb
        67      IF (dchisq.GT.0) THEN
        68          h1a= h1b
        69      ELSE IF (dchisq.LE.0) THEN
        70          deltah= (-deltah*2)/100
        71      END IF
        72      IF (chisqb.LE.6000) EXIT
        73  10 END DO
        74  WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
        75  !*****************************************************************
        76 
        77  END PROGRAM assignment
        78 
        79  Read more: https://www.physicsforums.com
     
    First question...are you compiling with Fortran 77 or Fortran 90? I think 90.

    I say this because you have typos in the program that may be doing something you did not intend. Recall that Fortran 90 operates on entire arrays when instructed to do so, i.e., no index is passed.

    Using numbered linsting above:

    line 44 calculates i-th index of Cs and only this index
    line 45 performs an array operation!...is this what you want?
    line 46 performs an array operation.
    line 47 takes a single index from c

    what is the point of performing arrays operations if ultimately you seem to just need a single index?

    deltah is set to 0.001 at the end of do-loop 20

    line 52 sets h1b to 0.001 and
    line 53 sets chisqb to 0

    this means that do-loop 30 does exactly the same thing on every iteration!

    in do-loop 30, you seem to correctly use parenthesis and the i index in your operations, yet, in
    line 58 you assign a single scalar to the entire 'c' array (since you did not specify an index)
    line 59 will pick up the intended value and assign to chisqb

    line 70 may flip the sign of deltah and divide it by 50, but this new value will be used only once at the beginning of do-loop 20; is this what you want? a short step back? just asking
     
  4. Jan 26, 2014 #3
    Hi! Thanks for the reply :) I'm compiling with a program called Force2.0, but it's definitely meant to be fortran 77. With skills like mine, wouldn't be surprised if I've screwed up a few times!

    I'll bullet point responses to keep things clear(ish).
    (44)- Riiiight. So, having i=1,17 there isn't running through all three equations and then summing? Crap.
    (45,46,47)- Yeah I think so- this whole do block is supposed to be calculating chi-squared for values of t(i)=1,17 (ie all seventeen values of t read in from the data file), then adding them all together to find the actual chi-squared for that particular value of lamda 1 (h1 as I've entered it).

    NB: I can't believe I didn't notice that the middle two equations are missing their array assignments. Dumbass. All four should be using arrays.

    - Yeah I wasn't sure of the placement of deltah at the time of posting, I've moved it outside of the main loop since then. I'll post my newest code at the end of this reply. Likewise with h1b- that's been moved out the way too. I've left chisqb there as it makes sense in my head that chisq should be zeroed each time.

    - Another delightful cockup on my part, thanks for pointing it out. These missing arrays must be from when I imported them from my original attempt (tried using functions and subroutines, did not end well), but still can't believe I (nor my friends who had a quick look) didn't notice them.

    - hmm no. The idea was it's supposed to run through the whole loop (only operates in the second inner loop in the new code), knock the step back and downsize, but only for that iteration (unless dchisq is negative again). I've edited it to take two steps back now and decrease step size further (though I may change the size of the decrease once the minimisation loop is working).

    Here's the newest code (I've added back in the missing array identifiers).

    Code (Text):
          PROGRAM assignment
         
          ! A program designed to fit experiemental data, using the method
          ! of least squares to minimise the associated chi-squared and
          ! obtain the four control parameters A,B,h1 and h2.
          !*****************************************************************

          IMPLICIT NONE
          INTEGER i,t0
          DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
          DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
          DOUBLE PRECISION deltah,Cs(17)

            OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
               DO i=1,17
                  READ(21,*)t(i),Ct(i),eCt(i)
               END DO
            CLOSE(21)

          !Read in data.txt as three one dimensional arrays.
          !*****************************************************************
          !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
          !DO i=1,17
          !   WRITE(21,*)t(i),Ct(i),eCt(i)
          !END DO
          !CLOSE(21)
          !
          !Just to check input file is being read correctly.
          !*****************************************************************
          !****************************Parameters***************************
         
          Aa= 0
          Ba= 0
          h1a= 0
          h2a= 0

          !*****************************************************************
          !**********************Minimising Lamda1 (h1)*********************
          deltah= 0.0001
          h1b= 0.001
          h1a= 0
               DO 10
          chisqa= 0
                  DO 20 i= 1,17
                     h1a= h1a
                     Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
                     Ctdiff(i)= Ct(i) - Cs(i)
                     c(i)= Ctdiff(i)**2/eCt(i)**2
                     chisqa= chisqa + c(i)
     20           END DO
          !Use initial h1 value of 0 to calculate start-point chisq.

          chisqb= 0
                  DO 30 i= 1,17
                     h1b= h1b + deltah
                     Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
                     Ctdiff(i)= Ct(i) - Cs(i)
                     c(i)= Ctdiff(i)**2/eCt(i)**2
                     chisqb= chisqb + c(i)
     30           END DO
          !First-step h1 used to find competing chisq for comparison.
            WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
            WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
          !Prints the two calculated chisq values to screen.
                     dchisq= chisqa - chisqb
                  IF (dchisq.GT.0) THEN
                     h1a= h1b
                  ELSE IF (dchisq.LE.0) THEN
                     deltah= deltah-((deltah*2)/100)
                  END IF
                 
               IF (chisqb.LE.6618.681) EXIT
     10        END DO
               WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
          !*****************************************************************
               
          END PROGRAM assignment
    Thankyou for the help :)

    EDIT: I imagine I'm going to have to use a function for the chi-squared calculation aren't I?
     
    Last edited: Jan 26, 2014
  5. Jan 26, 2014 #4

    Mark44

    Staff: Mentor

    I added some comments where I think you're going wrong. My comments look like this: !!! ***
    In one comment, you have three assignment statements where you are using array variables, but without indexes. This is wrong. Everywhere you have an array variable in a calculation, it should have an index. E.g., it should be Ct(i), not just Ct.

    In my other comments, you have a line of code that might be wrong. It looks like you commented it out, but I'm not sure what you did.

    Your variable names are terrible! For example, Aa,Ab,Ba,Bb,h1a,h1b. Can't you come up with some names that are self-explanatory? Also, why does your code meander all across the page? This makes it harder to see which statements are inside IF blocks and DO loops and such.
     
  6. Jan 26, 2014 #5
    Hi! I'll reply using bullet points for ease of reading. As always, thanks for taking the time to read through my mess :)
    Comments first!

    - RE: commented out variables- yeah I've commented them out whilst I work on the single parameter minimisation. They'll come back in one by one once it's working.
    - You caught my cockup too. No clue how I managed to go the last few hours without noticing. Guess I just didn't expect there to be issues with those blocks and ignored them- foolish.


    As for variable names, well, they made sense to me at the time but if you'd prefer something clearer I can alter them. RE: meandering code, I assume you mean the equation blocks? They weren't like that before, got whinged at by some here when they had a look so changed it, but I'll change it back.
    If you mean the code in general, I thought it was okay. Will alter it if you have any specific preferences.

    EDIT: Sorry about the thread thing.
     
    Last edited: Jan 26, 2014
  7. Jan 26, 2014 #6

    Mark44

    Staff: Mentor

    When you post code, it's good practice here at PF to put [code] ... [/code] tags around it. I've done this for you in your code.
    Probably be a good idea.
    Is your code producing the right results now?
    Don't do it for me - do it for yourself. When you review this code in a month or two, having meaningless variable names will make it harder for you to understand the code, even though you wrote it yourself.
    We call them assignments statements.
    Minor nit, but this Greek letter - λ - is written "lambda".
     
  8. Jan 26, 2014 #7
    Done!

    So as the assignment statements stand right now, are they not running through the calculation for every array labelled parameter using the i=1,17 statement? If no, it'll be better to move all of the assignment statements to a function and just execute that 17 times, yes? I'm not great with functions; going to throw one together quick and post it in an edit here, help with any needed corrections would be appreciated.

    No, still stuck at two lamda values.

    This is true. What might be a great idea at the time might not be later. Not particularly creative, but I've just renamed them after the do loop they're used in.
    I know re: lamda, but I didn't want "lamda1, lamda2" lengthening lines out. I usually use h in its place in fortran, since the actual character isn't available.
     
  9. Jan 26, 2014 #8
    I've written this function to replace the assignment statements in the two interior loops, leaving the summation statements in.

    Code (Text):

                  DO 20 i= 1,17
                     h1a= h1a
                     chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),Cs(i),eCt(i))
     20           END DO




          DOUBLE PRECISION FUNCTION c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)
          IMPLICIT NONE
          DOUBLE PRECISION Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt

                     Cs= exp(-h1loop1*t)!*Aa !+ Ba*exp(-h2loop1*t(i))
                     Ctdiff= Ct - Cs
                     c= Ctdiff**2/eCt**2

          RETURN
          END
    Will this do the trick?

    EDIT: It doesn't like the array indexes in the function call-in in the main program body. Not quite sure how to properly call it in...do I just assign the indexes inside the function itself (swear somebody told me you couldn't use them in functions) and write the call-in as c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)?
     
    Last edited: Jan 26, 2014
  10. Jan 27, 2014 #9

    Mark44

    Staff: Mentor

    There are several things wrong with what you're doing. The DO loop is presumably part of your main program. It would be in some code that looks like this:
    Code (Text):

    PROGRAM Whatever

    ! Declarations
       ...
       DO 20 i= 1,17
           h1a= h1a
           chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),Cs(i),eCt(i))
    20 END DO
       ...
    END PROGRAM Whatever
     
    In the code above you have a completely useless assignment statement: h1a = h1a. Here you are assigning the value of h1a to itself. That value is already there, so there's no point in doing this.

    In your main program you need a declaration for your c function. BTW, you really should have a better name for it than 'c'. For a function, a reasonable name is a noun that describes the value that is being returned.

    Below the main program are any subprograms, including your function, which is one type of subprogram (subroutine being the other).

    The main program has the statement that is the call to the function. (In 40+ years I've never heard anybody refer to it as a 'call-in'.) Your function call looks OK - it has the right number of parameters, and they appear to be the correct types. After you have your function declaration in place, I think things should work properly.

    It's lamBda - there's a 'B' in the name.
     
  11. Jan 27, 2014 #10
    I wasn't sure if h1a would need stating before the loop or not, now I do.
    c works for me. It's an intermediate value that has no use outside of its function, and it's only definition is as the sum of c(i) array- Chi-squared.
    I like to be different. That's what it does so works for me to call it that, I'm a physicist first not a programmer.

    Here's the code as it now stands.

    Code (Text):
          PROGRAM assignment

          ! A program designed to fit experiemental data, using the method
          ! of least squares to minimise the associated chi-squared and
          ! obtain the four control parameters A,B,h1 and h2.
          !*****************************************************************

          IMPLICIT NONE
          INTEGER i,t0
          DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17)
          DOUBLE PRECISION Aloop1,Aloop2,Bloop1,Bloop2,h1loop1,h1loop2
          DOUBLE PRECISION deltah,Cs(17),h2loop1,h2loop2
          DOUBLE PRECISION chisqa,chisqb,dchisq,c(17,17,17,17,17,17,17)

            OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
               DO i=1,17
                  READ(21,*)t(i),Ct(i),eCt(i)
               END DO
            CLOSE(21)

          !Read in data.txt as three one dimensional arrays.
          !*****************************************************************
          !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
          !DO i=1,17
          !   WRITE(21,*)t(i),Ct(i),eCt(i)
          !END DO
          !CLOSE(21)
          !
          !Just to check input file is being read correctly.
          !*****************************************************************
          !****************************Parameters***************************

          Aloop1= 0
          Bloop2= 0
          h1loop1= 0
          h2loop2= 0

          !*****************************************************************
          !**********************Minimising Lamda1 (h1)*********************
          deltah= 0.0001
          h1loop2= 0.001
          h1loop1= 0
               DO 10
          chisqa= 0
                  DO 20 i= 1,17
                     chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
         >           Cs(i),eCt(i))
     20           END DO
            WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1loop1
          !Use initial h1 value of 0 to calculate start-point chisq.
         
          chisqb= 0
                  DO 30 i= 1,17
                     h1loop2= h1loop2 + deltah
                     chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
         >           Cs(i),eCt(i))
     30           END DO
          !First-step h1 used to find competing chisq for comparison.
            WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1loop1
            WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1loop2
          !Prints the two calculated chisq values to screen.
                     dchisq= chisqa - chisqb
                  IF (dchisq.GT.0) THEN
                     h1loop1= h1loop2
                  ELSE IF (dchisq.LE.0) THEN
                     deltah= deltah-((deltah*2)/100)
                  END IF

               IF (chisqb.LE.6618.681) EXIT
     10        END DO
               WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1loop2
          !*****************************************************************

          END PROGRAM assignment

          !*****************************************************************
          !***************************FUNCTIONS*****************************
          DOUBLE PRECISION FUNCTION c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)
          IMPLICIT NONE
          DOUBLE PRECISION Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt

                     Cs= exp(-h1loop1*t)!*Aa !+ Ba*exp(-h2loop1*t(i))
                     Ctdiff= Ct - Cs
                     c= Ctdiff**2/eCt**2

          RETURN
          END FUNCTION c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)
    Just took a stab at defining c at the start of the program, wasn't entirely sure but there you go. This presents the following errors on compiling.


    test.f:46.36:

    chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:46.42:

    chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:46.50:

    chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    C:\Users\Nathan\Desktop\Documents\University of Birmingham MSc Physics and Technology of Nuclear Reactors\Lectures\fortran\test.f:46.55:

    chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:46.65:

    chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:46.71:

    chisqa= chisqa + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:47.23:

    > Cs(i),eCt(i))
    1
    Warning: Extension: REAL array index at (1)
    test.f:55.36:

    chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:55.42:

    chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:55.50:

    chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:55.55:

    chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:55.65:

    chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:55.71:

    chisqb= chisqb + c(Cs(i),h1loop1,t(i),Ctdiff(i),Ct(i),
    1
    Warning: Extension: REAL array index at (1)
    test.f:56.23:

    > Cs(i),eCt(i))
    1
    Warning: Extension: REAL array index at (1)
    test.f:1.72:

    PROGRAM assignment
    1
    test.f:78.72:

    DOUBLE PRECISION FUNCTION c(Cs,h1loop1,t,Ctdiff,Ct,Cs,eCt)
    2
    Error: Two main PROGRAMs at (1) and (2)
    test.f:13: error: size of variable 'c' is too large

    Surprised by the last one. I thought fortran 77 could handle seven dimensions to an array. Though perhaps I've just not entered something correctly.
     
  12. Jan 27, 2014 #11
    O.k., stay away from functions, for now. But please, just because you are getting help here, don't stop studying Fortran...you really need to for this assignment and as a physicist!

    It seems that your comments in the source apply to the code above them...that is not how it is typically done; instead, you first tell the reader what you are about to do, then you do it.

    From what I can tell, the only variables that need to be arrays are the ones you read from the file; so, I cleaned that.

    Deleted unnecessary stuff to force you to add it later with consistent variable naming.

    Changed the program a bit to force you to re-read it and focus on the correctness, once again.

    And that is it, I am not going to bother to explain EVERYTHING that is wrong with your latest code (you really went backwards); instead, please take a good look at the cleaned-up code that I post below...it at least compiles; but, as simple as it might be, I have not check the math at all or the looping or anything, I just looked at the Fortran.

    Code (Text):

          !
          ! A program designed to fit experiemental data, using the method
          ! of least squares to minimise the associated chi-squared and
          ! obtain the four control parameters A,B,h1 and h2.      
          !
          PROGRAM assignment
         
          IMPLICIT NONE
          INTEGER i
          DOUBLE PRECISION t(17), Ct(17), eCt(17)
          DOUBLE PRECISION h1loop1, h1loop2, deltah, Cs
          DOUBLE PRECISION chisqa, chisqb, dchisq

          OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
          DO i=1,17
              READ(21,*) t(i), Ct(i), eCt(i)
          END DO
          CLOSE(21)

          !**********************Minimising Lamda1 (h1)*********************
          deltah  = 0.0001
          h1loop2 = 0.001
          h1loop1 = 0.0 !Use initial value of 0 to calculate start-point chisq
          DO 10
              chisqa = 0.0
              DO 20 i = 1, 17
                  Cs     = exp(-h1loop1*t(i))
                  chisqa = chisqa + ( (Ct(i) - Cs)/eCt(i) )**2
     20       END DO
         
              chisqb = 0.0
              DO 30 i = 1, 17
                  h1loop2 = h1loop2 + deltah
                  Cs      = exp(-h1loop2*t(i))
                  chisqb  = chisqb + ( (Ct(i) - Cs)/eCt(i) )**2
     30       END DO

             !Print the two calculated chisq values to screen.
              WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1loop1
              WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1loop2
             
              dchisq = chisqa - chisqb
              IF (dchisq.GT.0.0) THEN
                  h1loop1 = h1loop2
              ELSE
                  deltah = deltah-((deltah*2)/100)
              END IF

              IF (chisqb.LE.6618.681) EXIT
     10   END DO
     
          WRITE(6,*) 'Chi-squared is', chisqb,' for Lamda1 = ', h1loop2
          END PROGRAM assignment
     
     
  13. Jan 27, 2014 #12
    lol oh don't worry, I'm going to be getting a lot more familiar with fortran77 over the next few years. We're only really getting into it now, but it'll be used in my lab reports later in the semester and my dissertation over the summer. Then, I hear reactor control is still done in fortran77...there'll be no escape!

    No probs, I'll move the comments around. I'll compare the new code to my own latest and learn from it, thanks. Updated my code, the iteration is still getting stuck between two values. This problem is ridiculous.
     
    Last edited: Jan 27, 2014
  14. Jan 27, 2014 #13

    SteamKing

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    You have defined 'c' as a DOUBLE PRECISION array in the main program and then again as a DOUBLE PRECISION function. Which is it? That's why you are getting a bunch of compiler errors.
     
  15. Jan 27, 2014 #14
    I would recommend reducing the number of points from 17 down to 3 and performing a hand calculation...I mean, really, hand calculate values one iteration at a time to really see what your equations are yielding, in particular, when to increment h1loop1 and h1loop2.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Fortran77 woes- single parameter least squares minimisation stuck!
  1. Python woes (Replies: 7)

Loading...