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
tresid_(n,m,d,e,matrixZ,mapZ,eval,iwork,work,res,info)38 void tresid_( n, m, d, e, matrixZ, mapZ, eval, iwork, work, res, info)
39 Integer *n, *m, *mapZ, *iwork, *info;
40 DoublePrecision *d, *e, *matrixZ, *eval, *work, *res;
41
42 /*
43
44 FORTRAN wrapper for tresid().
45
46 this subroutine computes the residual
47
48 res = max_{i} | T z_{i} - \lambda_{i} z_{i} |/( | T | * ulp )
49
50 where T is an n-by-n tridiagonal matrix,
51 ( \lambda_{i} , z_{i} ) is a standard eigen-pair, and
52 ULP = (machine precision) * (machine base)
53
54 |T| is the 1-norm of T
55 |T z_{i} .... | is the infinity-norm
56 res is reasonable if it is of order about 50 or less.
57
58 Arguments
59 ---------
60
61 n = (input) dimension of the matrix T
62 m = (input) number of eigenvalues/eigenvectors
63
64 d(1:n) = (input) diagonal of T
65
66 e(1) = junk
67 e(2:n) = (input) sub-diagonal of T
68 = super-diagonal of T
69
70 colZ = (input) DoublePrecision eigenvectors in packed storage
71 mapZ = (input) integer array of length m holding the id of the processors
72 holding the eigenvectors in colZ.
73 eval = (input ) array of m DoublePrecision eigenvalues
74
75 iwork = (workspace) array of Integer workspace
76 work = (workspace) array of DoublePrecision workspace
77
78 res = (output ) the residual described above.
79
80 info = (output ) = 0, then everything ok
81 < 0, then -info-th input variable is not valid
82
83
84 */
85 {
86
87 Integer ll, k, i;
88 Integer nvecsZ, me;
89
90 DoublePrecision *scratch;
91 DoublePrecision **buff_ptr, **ptrZ;
92
93
94 extern Integer mxmynd_();
95 extern Integer count_list();
96 extern void resid();
97
98 /*
99 usual story about error handling
100 */
101
102 ll = *n;
103
104 if ( ll < 1 )
105 return; /* error */
106
107 me = mxmynd_();
108
109 scratch = work;
110 buff_ptr = (DoublePrecision ** ) work;
111
112 nvecsZ = count_list( me, mapZ, m );
113
114 ptrZ = buff_ptr;
115 buff_ptr += nvecsZ + 1;
116 scratch += nvecsZ + 1;
117
118 k = 0;
119 for ( i = 0; i < nvecsZ; i++ ) {
120 ptrZ[i] = &matrixZ[k];
121 k += ll;
122 }
123
124 tresid( n, m, d, e, ptrZ, mapZ, eval, iwork, scratch, res, info);
125
126 return;
127 }
128
129