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
36 #include "globalp.c.h"
37
38 #define MSG_START 25000
39 #define ZERO ((DoublePrecision) 0.0e0)
40
41
42
residual_(n,matrixA,mapA,matrixB,mapB,m,matrixZ,mapZ,eval,iwork,work,res,info)43 void residual_( n, matrixA, mapA, matrixB, mapB, m, matrixZ, mapZ, eval, iwork, work, res, info)
44 Integer *n, *mapA, *mapB, *m, *mapZ, *iwork, *info;
45 DoublePrecision *matrixA, *matrixB, *matrixZ, *eval, *work, *res;
46
47 /*
48
49 this subroutine computes the residual
50
51 res = max_{i} | A z_{i} - \lambda_{i} B z_{i} |/( | A | * ulp )
52
53 where (lambda, z) is a generalized eigen-pair.
54 ULP = (machine precision) * (machine base)
55
56 |A z_{i} ... | is the infinity-norm,
57 |A| is the 1-norm of A,
58 res is reasonable if it is of order about 50 or less.
59
60
61 Assume that A and B are n x n symmetric matrices
62 in packed storage format column-distributed
63 across all the participating processors
64
65 on input:
66
67 n = (input) dimension of the matrix A, B
68 mapZ = (input) integer array of length m holding the processor
69 storing the matrices.
70 mapA, mapB = (input) integer array of length n holding the processor
71 storing the matrices.
72
73 matrixA, matrixB, matrixZ = (input) pointers to the appropriate locations
74 eval = (input, unchanged ) array of m DoublePrecision precision eigenvalues
75
76 m = number of columns in vecZ
77
78 work = scratch array
79
80 */
81 {
82
83 Integer ll, k, i, *iscrat, *mapvecA, *mapvecB, *mapvecZ;
84 Integer nvecsA, nvecsB, nvecsZ;
85 Integer me;
86
87 DoublePrecision *scratch;
88 DoublePrecision **buff_ptr, **ptrA, **ptrB, **ptrZ;
89
90
91 /*
92 blas call
93 */
94
95 extern Integer mxwrit_(), mxread_();
96 extern Integer menode_(), mxmynd_();
97
98 extern Integer fil_mapvec_();
99 extern Integer count_list();
100 extern void residual();
101
102 /*
103 usual story about error handling
104 */
105
106
107
108 ll = *n;
109
110 if ( ll < 1 )
111 return; /* error */
112
113 /*
114 should perform global sync to check for errors
115 */
116
117 me = mxmynd_();
118
119 scratch = work;
120 buff_ptr = (DoublePrecision ** ) work;
121
122 mapvecA = iwork;
123 iscrat = iwork;
124 nvecsA = fil_mapvec_( &me, &ll, mapA, mapvecA);
125 iscrat += nvecsA;
126
127 mapvecB = iscrat;
128 nvecsB = fil_mapvec_( &me, &ll, mapB, mapvecB);
129 iscrat += nvecsB;
130
131 mapvecZ = iscrat;
132 nvecsZ = fil_mapvec_( &me, m, mapZ, mapvecZ );
133 iscrat += nvecsZ;
134
135 ptrA = buff_ptr;
136 buff_ptr += nvecsA;
137 scratch += nvecsA + 1;
138
139 k = 0;
140 for ( i = 0; i < nvecsA; i++ ) {
141 ptrA[i] = &matrixA[k];
142 k += ( ll - mapvecA[i] );
143 }
144
145 ptrB = buff_ptr;
146 buff_ptr += nvecsB;
147 scratch += nvecsB + 1;
148
149 k = 0;
150 for ( i = 0; i < nvecsB; i++ ) {
151 ptrB[i] = &matrixB[k];
152 k += ( ll - mapvecB[i] );
153 }
154
155 ptrZ = buff_ptr;
156 buff_ptr += nvecsZ;
157 scratch += nvecsZ + 1;
158
159 k = 0;
160 for ( i = 0; i < nvecsZ; i++ ) {
161 ptrZ[i] = &matrixZ[k];
162 k += ll;
163 }
164
165 scratch += 3 * ( nvecsZ + 1 );
166
167 residual( n, ptrA, mapA, ptrB, mapB, m, ptrZ, mapZ,
168 eval, buff_ptr, iwork, scratch, res, info);
169
170 return;
171 }
172
173