/
ch2inv.f
74 lines (74 loc) · 2.26 KB
/
ch2inv.f
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
c-----------------------------------------------------------------------
c
c R : A Computer Langage for Statistical Data Analysis
c Copyright (C) 1996, 1997 Robert Gentleman and Ross Ihaka
c Copyright (C) 2001 The R Development Core Team
c
c This program is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 2 of the License, or
c (at your option) any later version.
c
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
c GNU General Public License for more details.
c
c You should have received a copy of the GNU General Public License
c along with this program; if not, write to the Free Software
c Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
c
c-----------------------------------------------------------------------
c
c ch2inv(x, ldx, n, v, info) computes the inverse of a positive-definite
c symmetric matrix from its choleski factorization.
c This can be used, e.g., to compute the dispersion matrix for the estimated
c parameters in a regression analysis.
c
c On Entry
c
c x double precision(ldx,n) {ldx >= n}
c the choleski decomposition or the qr decomposition
c as computed by dqrdc or dqrdc2. Only the
c UPPER TRIANGULAR part, x(i,j), 1 <= i <= j <= n, is accessed.
c
c ldx integer, ldx >= n
c the leading dimension of the array x
c
c n integer
c the number of rows and columns of the matrix x
c
c On Return
c
c v double precision(n,n)
c the value of inverse(x'x)
c
c This version dated Aug 24, 1996.
c Ross Ihaka, University of Auckland.
c
subroutine ch2inv(x, ldx, n, v, info)
c implicit none
integer n, ldx, info
double precision x(ldx, n), v(n, n)
c
double precision d
integer i, j, im1
c
do 20 i=1,n
if(x(i,i) .eq. 0.0d0) then
info = i
return
end if
do 10 j=i,n
v(i,j) = x(i,j)
10 continue
20 continue
call dpodi(v, n, n, d, 1)
do 40 i=1,n
im1 = i-1
do 30 j=1,im1
v(i,j) = v(j,i)
30 continue
40 continue
return
end