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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
C Copyright (c) 2003-2010 University of Florida
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 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 The GNU General Public License is included in this distribution
C in the file COPYRIGHT.
subroutine remove_diagonal(array_table,
* narray_table, index_table,
* nindex_table, segment_table, nsegment_table,
* block_map_table, nblock_map_table,
* scalar_table, nscalar_table,
* address_table, op)
c--------------------------------------------------------------------------
c
c The diagonal elements of an array(Array1) are removed and the array
c without the diagonal elements is put into another array(Array2).
c
c The instruction is executed as follows:
c
c execute remove_diagonal Array1 Array2
c
c where Array1 is the full input matrix, usually the fock matrix.
c Array2 is the matrix with zero diagonal elements.
c
c Both Array1 and Array2 must be "static" arrays.
c--------------------------------------------------------------------------
implicit none
include 'interpreter.h'
include 'trace.h'
include 'epsilon.h'
#ifdef ALTIX
include 'sheap.h'
#endif
integer narray_table, nindex_table, nsegment_table,
* nblock_map_table
integer op(loptable_entry)
integer array_table(larray_table_entry, narray_table)
integer index_table(lindex_table_entry, nindex_table)
integer segment_table(lsegment_table_entry, nsegment_table)
integer block_map_table(lblock_map_entry, nblock_map_table)
integer nscalar_table
double precision scalar_table(nscalar_table)
integer*8 address_table(narray_table)
integer*8 iarray, ievec, get_index_from_base
integer i, j, n
integer array, evec_array
integer array_type, evec_type
integer nindex, nindex_evec
integer ind(mx_array_index)
integer junk
double precision x(1)
#ifdef ALTIX
pointer (dptr, x)
#else
common x
#endif
if (dryrun) return
#ifdef ALTIX
dptr = dshptr
#endif
c----------------------------------------------------------------------------
c Locate the data for the input matrix.
c---------------------------------------------------------------------------
array = op(c_result_array)
evec_array = op(c_op1_array)
if (array .eq. 0 .or. evec_array .eq. 0) then
print *,'Error: remove_diagonal routine requires 2
* array arguments.'
print *,(op(i),i=1,loptable_entry)
call abort_job()
endif
array_type = array_table(c_array_type, array)
evec_type = array_table(c_array_type, evec_array)
if (array_type .ne. static_array .or.
* evec_type .ne. static_array) then
print *,'Error: Both arrays in remove_diagonal routine
* must be static.'
call abort_job()
endif
nindex = array_table(c_nindex, array)
nindex_evec = array_table(c_nindex, evec_array)
if (nindex .ne. nindex_evec) then
print *,'Error: Both arrays in remove_diagonal must use ',
* 'the same number of indices.'
print *,'First array has ',nindex,' indices.'
print *,'Second array has ',nindex_evec,' indices.'
call abort_job()
endif
do i = 1, nindex
ind(i) = array_table(c_index_array1+i-1,array)
enddo
c---------------------------------------------------------------------------
c Look up each array's address.
c---------------------------------------------------------------------------
iarray = get_index_from_base(address_table(array), x, 2)
ievec = get_index_from_base(address_table(evec_array), x, 2)
n = index_table(c_index_size, ind(1)) ! pick up length of index
if (nindex .eq. 2) call rdiag2(x(iarray),x(ievec),array,n)
return
end
subroutine rdiag2(array1,array2,array,n)
c---------------------------------------------------------------------------
c
c The diagonal elements of the array1 are removed and the output put into
c array2. The diagonal elements of the Fock matrix are also put into
c the epsilon array.
c
c---------------------------------------------------------------------------
implicit none
include 'epsilon.h'
integer n
double precision array1(n,n)
double precision array2(n,n)
integer a, b
integer array
do b = 1, n
do a = 1, n
if (a .ne. b) array2(a,b) = array1(a,b)
if (a .eq. b) array2(a,b) = 0.0
c write(6,*) a, b, array1(a,b), array2(a,b)
enddo
enddo
return
end
|
