1c-----------------------------------------------------------------------
2c
3c  R : A Computer Language for Statistical Data Analysis
4c  Copyright (C) 2003--2019  The R Foundation
5c  Copyright (C) 1996, 1997  Robert Gentleman and Ross Ihaka
6c
7c  This program is free software; you can redistribute it and/or modify
8c  it under the terms of the GNU General Public License as published by
9c  the Free Software Foundation; either version 2 of the License, or
10c  (at your option) any later version.
11c
12c  This program is distributed in the hope that it will be useful,
13c  but WITHOUT ANY WARRANTY; without even the implied warranty of
14c  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15c  GNU General Public License for more details.
16c
17c  You should have received a copy of the GNU General Public License
18c  along with this program; if not, a copy is available at
19c  https://www.R-project.org/Licenses/
20c
21c-----------------------------------------------------------------------
22c
23c   lminfl computes basic quantities useful for computing
24c   regression diagnostics.
25c
26c   on entry
27c
28c       x       double precision(ldx,k)
29c               the qr decomposition as computed by dqrdc or dqrdc2.
30c
31c       ldx     integer
32c               the leading dimension of the array x.
33c
34c       n       integer
35c               the number of rows of the matrix x.
36c
37c       k       integer
38c               the number of columns in the matrix x.
39c
40c       q       integer
41c               the number of columns of the (multivariate) y[]
42c
43c       qraux   double precision(k)
44c               auxiliary information about the qr decomposition.
45c
46c       resid   double precision(n, q)
47c               the residuals from the (possibly multivariate) regression.
48c
49c   on return
50c
51c       hat     double precision(n)
52c               the diagonal of the hat matrix.
53c
54c       sigma   double precision(n, q)
55c               the (i,j)-th element of sigma contains an estimate
56c               of the residual standard deviation for the y[,j]-model with
57c               the i-th case omitted.
58c
59c   Initial version dated Aug 24, 1996.
60c   Ross Ihaka, University of Auckland.
61c   'docoef' option added Feb 17, 2003;  Martin Maechler ETH Zurich.
62c   Handle hat == 1 case, Nov 2005.
63c   Argument 'tol' was real not double precision, Aug 2007
64c   'q' for multivariate (mlm) models added Sep, 2018;  Martin Maechler.
65c    docoef replaced by R code and removed Nov 2019.
66
67      subroutine lminfl(x, ldx, n, k, q, qraux, resid,
68     +     hat, sigma, tol)
69      integer ldx, n, k, q
70      double precision x(ldx,k), qraux(k), resid(n,q),
71     +     hat(n), sigma(n,q), tol
72
73      integer i, j, info
74      double precision sum, denom, dummy(1)
75c
76c     hat matrix diagonal h_ii = hat(i)    [sigma(i,1) as auxiliary] :
77c
78      do i = 1,n
79        hat(i) = 0.0d0
80      end do
81
82      do j = 1,k
83        do i = 1,n
84          sigma(i,1) = 0.0d0
85        end do
86        sigma(j,1) = 1.0d0
87        call dqrsl(x, ldx, n, k, qraux, sigma, sigma, dummy,
88     .             dummy, dummy, dummy, 10000, info)
89        do i = 1, n
90          hat(i) = hat(i) + sigma(i,1)*sigma(i,1)
91        end do
92      end do
93      do i = 1, n
94        if(hat(i) .ge. 1.0d0 - tol) hat(i) = 1.0d0
95      end do
96c
97c     estimated residual standard deviation  sigma(j,c)
98c
99      denom = (n - k - 1)
100      do j = 1,q
101         sum = 0.0d0
102         do i = 1,n
103            sum = sum + resid(i,j)*resid(i,j)
104         end do
105         do i = 1,n
106            if(hat(i) .lt. 1.0d0) then
107               sigma(i,j) = sqrt((sum -
108     +              resid(i,j)*resid(i,j)/(1.0d0-hat(i))) / denom)
109            else
110               sigma(i,j) = sqrt(sum/denom)
111            endif
112         end do
113      end do
114
115      return
116      end
117