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 <math.h>
36
37 #include "globalp.c.h"
38
39 /*
40 PeIGS internal routine, currently unsupported.
41
42 Let U = an upper triangular matrix
43 Let W = full matrix
44
45 This subroutine solves the matrix Y for the matrix problem
46
47 U Y = W
48
49 The matrix Y overwrites W.
50
51 Data structure:
52
53 The matrix U is row wrapped.
54 The matrix W is column wrapped.
55
56 */
57
upperUF_(n,rowU,mapU,nW,colW,mapW,iwork,work)58 void upperUF_ ( n, rowU, mapU, nW, colW, mapW, iwork, work)
59 Integer *n, *mapU, *nW, *mapW, *iwork;
60 DoublePrecision **rowU, **colW, *work;
61 {
62 /*
63 n = dimension of the matrix
64 mapU = index array holding the proces
65 rowU = DoublePrecision pointer to the array location of the i-th column
66
67 ditto for mapvecW, colW, mapW
68
69 iwork = integer scratch space of size p
70
71 remark: one can improve the x-fer by taking advantage size of the
72 vector by segmenting the column vectors
73
74 should put an info in the argument list with regard to underflow
75
76 */
77
78 static Integer IONE = 1;
79
80 Integer i, iii, k, me, isize;
81 Integer nvecsU, nvecsW;
82 Integer count_list(), i_U;
83 Integer nprocs;
84 Integer *mapvecU, *mapvecW;
85 Integer *iscrat;
86 Integer *proclist;
87
88 DoublePrecision *ptr;
89 DoublePrecision t, *buffer;
90 DoublePrecision dummy;
91 DoublePrecision safeulp;
92
93 extern void combine_vector();
94 extern Integer fil_mapvec_();
95 extern Integer indxL();
96
97 /*
98 blas calls
99 */
100
101 extern void dcopy_ ();
102 extern DoublePrecision ddot_ ();
103
104
105 /*
106 mxsubs communication calls
107 */
108
109 extern Integer mxwrit_ (), mxread_ (), mxmynd_ ();
110 extern void chol_pipe_bcast();
111 extern Integer mxmynd_ (), reduce_list2();
112
113 me = mxmynd_ ();
114
115 /*
116 info = 0;
117 */
118
119 safeulp = DLAMCHS;
120
121 iscrat = iwork;
122 mapvecU = iscrat;
123 nvecsU = fil_mapvec_ ( &me, n, mapU, mapvecU );
124
125 iscrat += nvecsU;
126 mapvecW = iscrat;
127 nvecsW = fil_mapvec_ ( &me, nW, mapW, mapvecW );
128 iscrat += nvecsW;
129
130 proclist = iscrat;
131 nprocs = reduce_list2( *nW, mapW, proclist);
132 iscrat += nprocs;
133
134 if ( nvecsU + nvecsW == 0 )
135 return;
136
137 buffer = work;
138
139 /*
140 count number of columns that I own
141 */
142
143
144 i_U = nvecsU-1;
145 for ( i = *n - 1 ; i >= 0; i-- ) {
146 isize = *n - i;
147 if ( mapU[i] == me ) {
148 /*
149 indx = indxL ( i, nvecsU, mapvecU );
150 */
151 dcopy_ ( &isize, rowU[i_U], &IONE, buffer, &IONE );
152 isize *= sizeof(DoublePrecision);
153 chol_pipe_bcast( (char *) buffer, isize, i, mapU[i], nprocs, proclist, iscrat );
154 i_U--;
155 }
156 else {
157 isize = *n - i;
158 isize *= sizeof(DoublePrecision);
159 iii = mapU[i];
160 chol_pipe_bcast( (char *) buffer, isize, i, iii, nprocs, proclist, iscrat);
161 }
162
163 /*
164 the buffer[i:n-1] now contains U[i,i:n-1]
165 */
166
167 *buffer = (DoublePrecision) 1.0 / *buffer;
168
169 /*
170 if ( *buffer > safeulp )
171 *buffer = (DoublePrecision) 1.0 / *buffer;
172 else {
173 *info = 1;
174 *buffer = 1.0e0/safeulp;
175 }
176
177 */
178 dummy = *buffer;
179 isize = *n-i-1;
180
181 for ( k = 0; k < nvecsW; k++ ) {
182 ptr = &colW[k][i + 1];
183 t = (DoublePrecision) -ddot_ ( &isize, &buffer[1], &IONE, ptr, &IONE);
184 t += colW[k][i];
185 t *= dummy;
186 colW[k][i] = t;
187 }
188 }
189
190 return;
191 }
192
193
194
195