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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
|
/* Factor.c */
#include "../Bridge.h"
#define MYDEBUG 1
#if MYDEBUG > 0
static int count_Factor = 0 ;
static double time_Factor = 0.0 ;
#endif
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------------------------
purpose -- to compute the factorization of A - sigma * B
note: all variables in the calling sequence are references
to allow call from fortran.
input parameters
data -- pointer to bridge data object
psigma -- shift for the matrix pencil
ppvttol -- pivot tolerance
*ppvttol = 0.0 --> no pivoting used
*ppvttol != 0.0 --> pivoting used, entries in factor are
bounded above by 1/pvttol in magnitude
output parameters
*pinertia -- on return contains the number of negative eigenvalues
*perror -- on return contains an error code
1 -- error found during factorization
0 -- normal return
-1 -- psigma is NULL
-2 -- ppvttol is NULL
-3 -- data is NULL
-4 -- pinertia is NULL
created -- 98aug10, cca & jcp
---------------------------------------------------------------------
*/
void
Factor (
double *psigma,
double *ppvttol,
void *data,
int *pinertia,
int *perror
) {
Bridge *bridge = (Bridge *) data ;
Chv *rootchv ;
ChvManager *chvmanager ;
double droptol=0.0, tau ;
double cpus[10] ;
int stats[20] ;
int nnegative, nzero, npositive, pivotingflag ;
#if MYDEBUG > 0
double t1, t2 ;
MARKTIME(t1) ;
count_Factor++ ;
fprintf(stdout, "\n (%d) Factor()", count_Factor) ;
fflush(stdout) ;
#endif
/*
---------------
check the input
---------------
*/
if ( psigma == NULL ) {
fprintf(stderr, "\n error in Factor()"
"\n psigma is NULL\n") ;
*perror = -1 ; return ;
}
if ( ppvttol == NULL ) {
fprintf(stderr, "\n error in Factor()"
"\n ppvttol is NULL\n") ;
*perror = -2 ; return ;
}
if ( data == NULL ) {
fprintf(stderr, "\n error in Factor()"
"\n data is NULL\n") ;
*perror = -3 ; return ;
}
if ( pinertia == NULL ) {
fprintf(stderr, "\n error in Factor()"
"\n pinertia is NULL\n") ;
*perror = -4 ; return ;
}
if ( perror == NULL ) {
fprintf(stderr, "\n error in Factor()"
"\n perror is NULL\n") ;
return ;
}
/*
----------------------------------
set the shift in the pencil object
----------------------------------
*/
bridge->pencil->sigma[0] = -(*psigma) ;
bridge->pencil->sigma[1] = 0.0 ;
/*
-----------------------------------------------------
clear the front matrix and submatrix mananger objects
-----------------------------------------------------
*/
FrontMtx_clearData(bridge->frontmtx);
SubMtxManager_clearData(bridge->mtxmanager);
/*
-----------------------------------------------------------
set the pivot tolerance.
NOTE: spooles's "tau" parameter is a bound on the magnitude
of the factor entries, and is the recipricol of that of the
pivot tolerance of the lanczos code
-----------------------------------------------------------
*/
if ( *ppvttol == 0.0 ) {
tau = 10.0 ;
pivotingflag = SPOOLES_NO_PIVOTING ;
} else {
tau = (1.0)/(*ppvttol) ;
pivotingflag = SPOOLES_PIVOTING ;
}
/*
----------------------------------
initialize the front matrix object
----------------------------------
*/
FrontMtx_init(bridge->frontmtx, bridge->frontETree, bridge->symbfacIVL,
SPOOLES_REAL, SPOOLES_SYMMETRIC, FRONTMTX_DENSE_FRONTS,
pivotingflag, NO_LOCK, 0, NULL, bridge->mtxmanager,
bridge->msglvl, bridge->msgFile) ;
/*
-------------------------
compute the factorization
-------------------------
*/
chvmanager = ChvManager_new() ;
ChvManager_init(chvmanager, NO_LOCK, 1);
IVfill(20, stats, 0) ;
DVfill(10, cpus, 0.0) ;
rootchv = FrontMtx_factorPencil(bridge->frontmtx, bridge->pencil, tau,
droptol, chvmanager, perror, cpus,
stats, bridge->msglvl, bridge->msgFile);
ChvManager_free(chvmanager);
/*
----------------------------
if matrix is singular then
set error flag and return
----------------------------
*/
if ( rootchv != NULL ) {
*perror = 1 ;
return ;
}
/*
------------------------------------------------------------------
post-process the factor matrix, convert from fronts to submatrices
------------------------------------------------------------------
*/
FrontMtx_postProcess(bridge->frontmtx, bridge->msglvl, bridge->msgFile);
/*
-------------------
compute the inertia
-------------------
*/
FrontMtx_inertia(bridge->frontmtx, &nnegative, &nzero, &npositive) ;
*pinertia = nnegative;
/*
------------------------------------------------------------------
set the error. (this is simple since when the spooles codes detect
a fatal error, they print out a message to stderr and exit.)
------------------------------------------------------------------
*/
*perror = 0 ;
#if MYDEBUG > 0
MARKTIME(t2) ;
time_Factor += t2 - t1 ;
fprintf(stdout, ", %8.3f seconds, %8.3f total time",
t2 - t1, time_Factor) ;
fflush(stdout) ;
#endif
return ; }
/*--------------------------------------------------------------------*/
|
