File: thybrj1.c

package info (click to toggle)
minpack 19961126-13
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 2,676 kB
  • ctags: 643
  • sloc: sh: 8,051; fortran: 2,400; ansic: 736; makefile: 137; awk: 13
file content (82 lines) | stat: -rw-r--r-- 1,839 bytes parent folder | download | duplicates (7)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
/*      driver for hybrj1 example. */

#include <stdio.h>
#include <math.h>
#include <minpack.h>

void fcn(int *n, double *x, double *fvec, double *fjac, int *ldfjac, 
	 int *iflag);

int main()
{
  int j, n, ldfjac, info, lwa;
  double tol, fnorm;
  double x[9], fvec[9], fjac[9*9], wa[99];
  int one=1;

  n = 9;

/*      the following starting values provide a rough solution. */

  for (j=1; j<=9; j++)
    {
      x[j-1] = -1.;
    }

  ldfjac = 9;
  lwa = 99;

/*      set tol to the square root of the machine precision. */
/*      unless high solutions are required, */
/*      this is the recommended setting. */

  tol = sqrt(dpmpar_(&one));

  hybrj1_(&fcn, &n, x, fvec, fjac, &ldfjac, &tol, &info, wa, &lwa);

  fnorm = enorm_(&n, fvec);

  printf("      final l2 norm of the residuals%15.7g\n\n", fnorm);
  printf("      exit parameter                %10i\n\n", info);
  printf("      final approximate solution\n");
  for (j=1; j<=n; j++) printf("%s%15.7g", j%3==1?"\n     ":"", x[j-1]);
  printf("\n");

  return 0;
}

void fcn(int *n, double *x, double *fvec, double *fjac, int *ldfjac, 
	 int *iflag)
{
  /*      subroutine fcn for hybrj1 example. */

  int j, k;
  double one=1, temp, temp1, temp2, three=3, two=2, zero=0, four=4;

  if (*iflag != 2)
    {
      for (k = 1; k <= *n; k++)
	{
	  temp = (three - two*x[k-1])*x[k-1];
	  temp1 = zero;
	  if (k != 1) temp1 = x[k-1-1];
	  temp2 = zero;
	  if (k != *n) temp2 = x[k+1-1];
	  fvec[k-1] = temp - temp1 - two*temp2 + one;
	}
    }
  else
    {
     for (k = 1; k <= *n; k++)
       {
	 for (j = 1; j <= *n; j++)
	   {
	     fjac[k-1 + *ldfjac*(j-1)] = zero;
	   }
         fjac[k-1 + *ldfjac*(k-1)] = three - four*x[k-1];
         if (k != 1) fjac[k-1 + *ldfjac*(k-1-1)] = -one;
         if (k != *n) fjac[k-1 + *ldfjac*(k+1-1)] = -two;
       }
    }
  return;
}