Skip to content

Commit 842bc0e

Browse files
HugoMValeHugoMVale
committed
fix scale_mat shape input
1 parent b3c4d06 commit 842bc0e

File tree

2 files changed

+138
-16
lines changed

2 files changed

+138
-16
lines changed

‎report5.txt‎

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,118 @@
1+
******************************************************************************
2+
* ODRPACK REPORT *
3+
******************************************************************************
4+
5+
6+
7+
*** Derivative checking report for fit by method of ODR ***
8+
9+
10+
For response 1 of observation 1
11+
12+
User
13+
Supplied Relative Derivative
14+
Derivative WRT Value Difference Assessment
15+
16+
BETA( 1) 1.63E+00 3.11E-08 Verified
17+
BETA( 2) 3.21E+00 1.75E-08 Verified
18+
DELTA( 1, 1) 1.63E+00 7.31E-09 Verified
19+
20+
Number of reliable digits in function results 16
21+
(estimated by ODRPACK)
22+
23+
Number of digits of agreement required between
24+
user supplied and finite difference derivative for
25+
user supplied derivative to be considered verified 4
26+
27+
Row number at which derivatives were checked 1
28+
29+
-values of the explanatory variables at this row
30+
31+
X( 1, 1) 9.82000000E-01
32+
******************************************************************************
33+
* ODRPACK REPORT *
34+
******************************************************************************
35+
36+
37+
*** Initial summary for fit by method of ODR ***
38+
39+
--- Problem Size:
40+
N = 4 (number with nonzero weight = 4)
41+
NQ = 1
42+
M = 1
43+
NP = 2 (number unfixed = 2)
44+
45+
--- Control Values:
46+
JOB = 01020
47+
= ABCDE, where
48+
A=0 ==> fit is not a restart.
49+
B=1 ==> deltas are initialized by user.
50+
C=0 ==> covariance matrix will be computed using
51+
derivatives re-evaluated at the solution.
52+
D=2 ==> derivatives are supplied by user.
53+
derivatives were checked.
54+
results appear correct.
55+
E=0 ==> method is explicit ODR.
56+
NDIGIT = 16 (estimated by ODRPACK)
57+
TAUFAC = 1.00E+00
58+
59+
--- Stopping Criteria:
60+
SSTOL = 1.49E-08 (sum of squares stopping tolerance)
61+
PARTOL = 3.67E-11 (parameter stopping tolerance)
62+
MAXIT = 45 (maximum number of iterations)
63+
64+
--- Initial Weighted Sum of Squares = 1.46854548E+05
65+
Sum of Squared Weighted Deltas = 0.00000000E+00
66+
Sum of Squared Weighted Epsilons = 1.46854548E+05
67+
68+
--- Function Parameter Summary:
69+
70+
Index BETA(K) Fixed Scale LOWER(K) UPPER(K)
71+
72+
(K) (IFIXB) (SCLB)
73+
74+
1 2.00E+00 NO 5.00E-01 0.00E+000 1.00E+001
75+
2 5.00E-01 NO 5.00E-01 0.00E+000 9.00E-001
76+
77+
--- Explanatory Variable and Delta Weight Summary:
78+
79+
Index X(I,J) DELTA(I,J) Fixed Scale Weight
80+
81+
(I,J) (IFIXX) (SCLD) (WD)
82+
83+
1,1 9.820E-01 0.000E+00 NO 1.66E-01 1.00E+00
84+
N,1 6.010E+00 0.000E+00 NO 1.66E-01 1.00E+00
85+
86+
--- Response Variable and Epsilon Error Weight Summary:
87+
88+
Index Y(I,L) Weight
89+
(I,L) (WE)
90+
91+
1,1 2.700E+00 1.000E+00
92+
N,1 4.030E+02 1.000E+00
93+
94+
*** Final summary for fit by method of ODR ***
95+
96+
--- Stopping Conditions:
97+
INFO = 1 ==> sum of squares convergence.
98+
NITER = 17 (number of iterations)
99+
NFEV = 64 (number of function evaluations)
100+
NJEV = 18 (number of jacobian evaluations)
101+
IRANK = 0 (rank deficiency)
102+
RCOND = 7.68E-02 (inverse condition number)
103+
ISTOP = 0 (returned by user from subroutine FCN)
104+
105+
--- Final Weighted Sums of Squares = 2.67371503E-01
106+
Sum of Squared Weighted Deltas = 2.46886137E-01
107+
Sum of Squared Weighted Epsilons = 2.04853651E-02
108+
109+
--- Residual Standard Deviation = 3.65630621E-01
110+
Degrees of Freedom = 2
111+
112+
--- Estimated BETA(J), J = 1, ..., NP:
113+
114+
BETA LOWER UPPER S.D. ___ 95% Confidence ___
115+
BETA Interval
116+
117+
1 1.63337602E+00 0.00E+00 1.00E+01 5.02E-01 -5.26E-01 to 3.79E+00
118+
2 9.00000000E-01 0.00E+00 9.00E-01 7.44E-02 5.80E-01 to 1.22E+00

‎src/odrpack_core.f90‎

Lines changed: 20 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1105,7 +1105,7 @@ subroutine eval_jac &
11051105
if (ifixb(1) < 0) then
11061106
do j = 1, np
11071107
call scale_mat(n, nq, we1, ldwe, ld2we, &
1108-
reshape(fjacb(:, j, :), [n, nq]), & ! this reshape should not be required. probably compiler bug
1108+
reshape(fjacb(:, j, :), [n, nq]), &
11091109
tempret(1:n, 1:nq))
11101110
fjacb(:, j, :) = tempret(1:n, 1:nq)
11111111
end do
@@ -1114,7 +1114,9 @@ subroutine eval_jac &
11141114
do j = 1, np
11151115
if (ifixb(j) >= 1) then
11161116
j1 = j1 + 1
1117-
call scale_mat(n, nq, we1, ldwe, ld2we, fjacb(:, j, :), tempret(1:n, 1:nq))
1117+
call scale_mat(n, nq, we1, ldwe, ld2we, &
1118+
reshape(fjacb(:, j, :), [n, nq]), &
1119+
tempret(1:n, 1:nq))
11181120
fjacb(:, j1, :) = tempret(1:n, 1:nq)
11191121
end if
11201122
end do
@@ -1123,7 +1125,9 @@ subroutine eval_jac &
11231125
! Weight the Jacobian's wrt DELTA as appropriate
11241126
if (isodr) then
11251127
do j = 1, m
1126-
call scale_mat(n, nq, we1, ldwe, ld2we, fjacd(:, j, :), tempret(1:n, 1:nq))
1128+
call scale_mat(n, nq, we1, ldwe, ld2we, &
1129+
reshape(fjacd(:, j, :), [n, nq]), &
1130+
tempret(1:n, 1:nq))
11271131
fjacd(:, j, :) = tempret(1:n, 1:nq)
11281132
end do
11291133
end if
@@ -5369,17 +5373,17 @@ pure subroutine scale_mat(n, m, wt, ldwt, ld2wt, t, wtt)
53695373
integer, intent(in) :: m
53705374
!! Number of columns of data in `t`.
53715375
real(wp), intent(in), target :: wt(..)
5372-
!! Array of shape `(ldwt,ld2wt,m)` holding the weights.
5376+
!! Array of shape conformable to `(ldwt,ld2wt,m)` holding the weights.
53735377
integer, intent(in) :: ldwt
53745378
!! Leading dimension of array `wt`.
53755379
integer, intent(in) :: ld2wt
53765380
!! Second dimension of array `wt`.
53775381
real(wp), intent(in), target :: t(..)
5378-
!! Array of shape `(n,m)` being scaled by `wt`.
5382+
!! Array of shape conformable to `(n,m)` being scaled by `wt`.
53795383
real(wp), intent(out), target :: wtt(..)
5380-
!! Array of shape `(n,m)` holding the result of weighting array `t` by `wt`. Array
5381-
!! `wtt` can be the same as `t` only if the arrays in `wt` are upper triangular with
5382-
!! zeros below the diagonal.
5384+
!! Array of shape conformable to `(n,m)` holding the result of weighting array `t` by
5385+
!! array `wt`. Array `wtt` can be the same as `t` only if the arrays in `wt` are upper
5386+
!! triangular with zeros below the diagonal.
53835387

53845388
! Local scalars
53855389
integer :: i, j
@@ -5391,6 +5395,13 @@ pure subroutine scale_mat(n, m, wt, ldwt, ld2wt, t, wtt)
53915395

53925396
if (n == 0 .or. m == 0) return
53935397

5398+
select rank (wt)
5399+
rank (1); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
5400+
rank (2); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
5401+
rank (3); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
5402+
rank default; error stop "Invalid rank of `wt`."
5403+
end select
5404+
53945405
select rank (t)
53955406
rank (1); t_(1:n, 1:m) => t
53965407
rank (2); t_(1:n, 1:m) => t
@@ -5400,14 +5411,7 @@ pure subroutine scale_mat(n, m, wt, ldwt, ld2wt, t, wtt)
54005411
select rank (wtt)
54015412
rank (1); wtt_(1:n, 1:m) => wtt
54025413
rank (2); wtt_(1:n, 1:m) => wtt
5403-
rank default; error stop "Invalid rank of `wt`."
5404-
end select
5405-
5406-
select rank (wt)
5407-
rank (1); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
5408-
rank (2); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
5409-
rank (3); wt_(1:ldwt, 1:ld2wt, 1:m) => wt
5410-
rank default; error stop "Invalid rank of `wt`."
5414+
rank default; error stop "Invalid rank of `wtt`."
54115415
end select
54125416

54135417
if (wt_(1, 1, 1) >= zero) then

0 commit comments

Comments
 (0)