1 /*
2 $Id$
3 *======================================================================
4 *
5 * DISCLAIMER
6 *
7 * This material was prepared as an account of work sponsored by an
8 * agency of the United States Government. Neither the United States
9 * Government nor the United States Department of Energy, nor Battelle,
10 * nor any of their employees, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
11 * ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY,
12 * COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT,
13 * SOFTWARE, OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT
14 * INFRINGE PRIVATELY OWNED RIGHTS.
15 *
16 * ACKNOWLEDGMENT
17 *
18 * This software and its documentation were produced with Government
19 * support under Contract Number DE-AC06-76RLO-1830 awarded by the United
20 * States Department of Energy. The Government retains a paid-up
21 * non-exclusive, irrevocable worldwide license to reproduce, prepare
22 * derivative works, perform publicly and display publicly by or for the
23 * Government, including the right to distribute to other Government
24 * contractors.
25 *
26 *======================================================================
27 *
28 * -- PEIGS routine (version 2.1) --
29 * Pacific Northwest Laboratory
30 * July 28, 1995
31 *
32 *======================================================================
33 */
34 #include <stdio.h>
35 #include <stdlib.h>
36
37 #include "globalp.c.h"
38
ortho(n,m,colZ,mapZ,ibuffptr,iwork,work,ort,info)39 void ortho( n, m, colZ, mapZ, ibuffptr, iwork, work, ort, info)
40 Integer *n, *m, *mapZ, *iwork, *info;
41 DoublePrecision **colZ, **ibuffptr, *work, *ort;
42
43 /*
44
45 This subroutine computes the infinity-norm measure:
46
47 ort = max_i | (Z^t.Z)_i - I_i | / ULP,
48
49 for the standard symmetric eigensystem problem where
50
51 Z is N-by-M
52 I is M-by-M.
53 Z_i is an eigenvector,
54 (Z^t.Z)_i is the i-th column of Z^t.Z
55 I_i is the i-th column of the m-by-m identity matrix
56 ULP = (machine precision) * (machine base)
57 |.| is the infinity-norm.
58
59 res is reasonable if it is of order about 50 or less.
60
61
62 MUST have M <= N. If M > N then program exits.
63 This is not a limitation of this subroutine as M > N
64 implies the columns of Z are linearly dependent, which
65 implies "ort" will always be large in this case.
66
67 Arguments
68 ---------
69 In the following:
70
71 INTEGER = "pointer to Integer"
72 DOUBLE PRECISION = "pointer to DoublePrecision"
73
74 me = this processor's id (= mxmynd_())
75 nprocs = number of allocated processors ( = mxnprc_())
76 nvecsZ = number of entries in mapZ equal to me
77 (= count_list( me, mapZ, m ))
78 nvecsZ_max = maximum number of entries in mapZ equal to i,
79 i = 0 to nprocs-1
80 (= max over i of count_list( i, mapZ, m ))
81 sDP = sizeof( DoublePrecision )
82
83
84 n....... (input) INTEGER
85 Number of rows in Z
86
87 m....... (input) INTEGER
88 number of columns in Z (i.e., # of
89 eigenvalues/eigenvectors).
90 Must have m <= n.
91
92 colZ ... (input) array of pointers to DoublePrecision,
93 length (nvecsZ)
94 The part of matrix Z owned by this processer, stored
95 in packed format, i.e., colZ[i] points to the start
96 of the i-th column of matrix Z owned by this
97 processor, i = 0 to nvecsZ-1.
98
99 mapZ ... (input) INTEGER array, length (m)
100 The i-th column of matrix Z is owned by processor
101 mapZ[i], i = 0 to m-1.
102
103 ibuffptr (workspace) array of pointers to DoublePrecision,
104 length(nvecsZ)
105
106 iwork... (workspace) INTEGER array, length( 7 * m )
107
108 work.... (workspace) DOUBLE PRECISION array,
109 length( nvecsZ * n + maximum( d_mxm25,
110 mxlbuf_() / sDP + 1 )
111 where d_mxm25 = maximum ( 2*m, (nvecsZ + 2*nvecsZ_max)*n )
112
113 ort..... (output) INTEGER
114 the residual described above.
115
116 info.... (output) INTEGER
117 = 0, not currently used
118
119 */
120 {
121
122 static Integer ONE = 1;
123
124 Integer ll, i, j, *iscrat, *mapvecZ;
125 Integer nvecsZ;
126 Integer me;
127 Integer nprocs, *proclist;
128
129 DoublePrecision ulp;
130 DoublePrecision *ptr, *scrat;
131 DoublePrecision **vecZ1; /* copies of the vecZ matrix */
132
133 /*
134 blas call
135 */
136
137 extern void dcopy_();
138 extern void mxm25();
139
140 extern Integer mxmynd_();
141
142 extern Integer fil_mapvec_();
143 extern Integer reduce_list2();
144 extern void mat_max();
145
146 /*
147 usual story about error handling
148 */
149
150 ll = *n;
151
152 *ort = 0.e0;
153 *info = 0;
154
155 if ( ll < 1 )
156 return; /* error */
157
158 me = mxmynd_();
159
160 if( *n < *m ) {
161 fprintf(stderr,
162 "Error in routine ortho m (=%d) < n (=%d). me = %d \n",
163 *m, *n, me);
164 exit(-1);
165 }
166
167 iscrat = iwork;
168 mapvecZ = iscrat;
169 nvecsZ = fil_mapvec_( &me, m, mapZ, mapvecZ );
170 iscrat += nvecsZ;
171
172 proclist = iscrat;
173 nprocs = reduce_list2( *m, mapZ, proclist);
174 iscrat += nprocs;
175
176 if ( nvecsZ == 0 )
177 return;
178
179 scrat = work;
180
181 vecZ1 = ibuffptr;
182
183 /*
184 copy over the matrix to a buffer area
185 */
186
187 ptr = scrat;
188 for ( i = 0; i < nvecsZ; i++ ) {
189 vecZ1[i] = ptr;
190 dcopy_( &ll, colZ[i], &ONE, ptr, &ONE );
191 ptr += ll;
192 }
193
194 scrat = ptr;
195
196
197 /*
198 Z1 is a copy of Z
199
200 Z1 <- Z^t . Z1;
201 */
202
203 mxm25 ( m, &ll, colZ, mapZ, m, vecZ1, mapZ, vecZ1, iscrat, scrat);
204
205 for ( i = 0; i < nvecsZ; i++ ) {
206 j = mapvecZ[i];
207 vecZ1[i][j] -= 1.0e0;
208 }
209
210 mat_max ( m, m, vecZ1, mapZ, ort, iscrat, scrat);
211
212 ulp = DLAMCHE * DLAMCHB;
213
214 *ort = *ort / ulp;
215
216 return;
217 }
218