1 /*
2 ARPACK++ v1.2 2/20/2000
3 c++ interface to ARPACK code.
4
5 MODULE ARCSPen.h.
6 Arpack++ class ARchSymMPencil definition.
7 (CHOLMOD wrapper)
8
9 Author of this class:
10 Martin Reuter
11 Date 11/05/2012
12
13 Arpack++ Author:
14 Francisco Gomes
15
16 ARPACK Authors
17 Richard Lehoucq
18 Danny Sorensen
19 Chao Yang
20 Dept. of Computational & Applied Mathematics
21 Rice University
22 Houston, Texas
23 */
24
25 #ifndef ARCSPEN_H
26 #define ARCSPEN_H
27
28 //#include "arch.h"
29 //#include "arerror.h"
30 #include "blas1c.h"
31 //#include "lapackc.h"
32 #include "arcsmat.h"
33
34
35 template<class ARTYPE>
36 class ARchSymPencil
37 {
38
39 protected:
40
41 ARchSymMatrix<ARTYPE>* A;
42 ARchSymMatrix<ARTYPE>* B;
43 cholmod_factor *LAsB ;
44 bool factoredAsB;
45 cholmod_common c ;
46
47 virtual void Copy(const ARchSymPencil& other);
48
49 // void SparseSaxpy(ARTYPE a, ARTYPE x[], int xind[], int nx, ARTYPE y[],
50 // int yind[], int ny, ARTYPE z[], int zind[], int& nz);
51
52 // void ExpandAsB();
53
54 // void SubtractAsB(ARTYPE sigma);
55
56 public:
57
IsFactored()58 bool IsFactored() { return factoredAsB; }
59
60 void FactorAsB(ARTYPE sigma);
61
MultAv(ARTYPE * v,ARTYPE * w)62 void MultAv(ARTYPE* v, ARTYPE* w) { A->MultMv(v,w); }
63
MultBv(ARTYPE * v,ARTYPE * w)64 void MultBv(ARTYPE* v, ARTYPE* w) { B->MultMv(v,w); }
65
66 void MultInvBAv(ARTYPE* v, ARTYPE* w);
67
68 void MultInvAsBv(ARTYPE* v, ARTYPE* w);
69
70 void DefineMatrices(ARchSymMatrix<ARTYPE>& Ap, ARchSymMatrix<ARTYPE>& Bp);
71
ARchSymPencil()72 ARchSymPencil() { factoredAsB = false; A=NULL; B=NULL; LAsB=NULL; cholmod_start (&c) ; }
73 // Short constructor that does nothing.
74
75 ARchSymPencil(ARchSymMatrix<ARTYPE>& Ap, ARchSymMatrix<ARTYPE>& Bp);
76 // Long constructor.
77
ARchSymPencil(const ARchSymPencil & other)78 ARchSymPencil(const ARchSymPencil& other) { cholmod_start (&c) ; Copy(other); }
79 // Copy constructor.
80
~ARchSymPencil()81 virtual ~ARchSymPencil() { if (LAsB) cholmod_free_factor(&LAsB,&c); cholmod_finish (&c) ;}
82 // Destructor.
83
84 ARchSymPencil& operator=(const ARchSymPencil& other);
85 // Assignment operator.
86
87 };
88
89 // ------------------------------------------------------------------------ //
90 // ARchSymPencil member functions definition. //
91 // ------------------------------------------------------------------------ //
92
93
94 template<class ARTYPE>
Copy(const ARchSymPencil<ARTYPE> & other)95 inline void ARchSymPencil<ARTYPE>::Copy(const ARchSymPencil<ARTYPE>& other)
96 {
97 if (LAsB) cholmod_free_factor(&LAsB,&c);
98 A = other.A;
99 B = other.B;
100 factoredAsB = other.factoredAsB;
101 if (factoredAsB)
102 LAsB = cholmod_copy_factor(other.LAsB,&c);
103
104 } // Copy.
105
106 /*
107 template<class ARTYPE>
108 void ARchSymPencil<ARTYPE>::
109 SparseSaxpy(ARTYPE a, ARTYPE x[], int xind[], int nx, ARTYPE y[],
110 int yind[], int ny, ARTYPE z[], int zind[], int& nz)
111 // A strongly sequential (and inefficient) sparse saxpy algorithm.
112 {
113
114 int ix, iy;
115
116 nz = 0;
117 if ((nx == 0) || (a == (ARTYPE)0)) {
118 copy(ny,y,1,z,1);
119 for (iy=0; iy!=ny; iy++) zind[iy] = yind[iy];
120 nz = ny;
121 return;
122 }
123 if (ny == 0) {
124 copy(nx,x,1,z,1);
125 scal(nx,a,z,1);
126 for (ix=0; ix!=nx; ix++) zind[ix] = xind[ix];
127 nz = nx;
128 return;
129 }
130 ix = 0;
131 iy = 0;
132 while (true) {
133 if (xind[ix] == yind[iy]) {
134 zind[nz] = xind[ix];
135 z[nz++] = a*x[ix++]+y[iy++];
136 if ((ix == nx)||(iy == ny)) break;
137 }
138 else if (xind[ix] < yind[iy]) {
139 zind[nz] = xind[ix];
140 z[nz++] = a*x[ix++];
141 if (ix == nx) break;
142 }
143 else {
144 zind[nz] = yind[iy];
145 z[nz++] = y[iy++];
146 if (iy == ny) break;
147 }
148 }
149 while (iy < ny) {
150 zind[nz] = yind[iy];
151 z[nz++] = y[iy++];
152 }
153 while (ix < nx) {
154 zind[nz] = xind[ix];
155 z[nz++] = x[ix++];
156 }
157
158 } // SparseSaxpy.
159
160
161 template<class ARTYPE>
162 void ARchSymPencil<ARTYPE>::ExpandAsB()
163 {
164
165 int i, j, k, n;
166 int *pcol, *irow, *index, *pos;
167 ARTYPE *value;
168
169 // Initializing variables.
170
171 n = AsB.n;
172 index = AsB.index;
173 value = AsB.value;
174 irow = &index[n+1];
175 pcol = new int[AsB.n+1];
176 pos = new int[AsB.n+1];
177 for (i=0; i<=n; i++) pcol[i] = index[i];
178 for (i=0; i<=n; i++) pos[i] = 0;
179
180 // Counting the elements in each column of AsB.
181
182 if (AsB.uplo == 'U') {
183
184 for (i=0; i!=n; i++) {
185 k = pcol[i+1];
186 if ((k!=pcol[i])&&(irow[k-1]==i)) k--;
187 for (j=pcol[i]; j<k; j++) pos[irow[j]]++;
188 }
189
190 }
191 else { // uplo == 'L'
192
193 for (i=0; i!=n; i++) {
194 k = pcol[i];
195 if ((k!=pcol[i+1])&&(irow[k]==i)) k++;
196 for (j=k; j<pcol[i+1]; j++) pos[irow[j]]++;
197 }
198
199 }
200
201 // Summing up index elements.
202
203 for (i=0; i<n; i++) pos[i+1] += pos[i];
204 for (i=n; i>0; i--) index[i] += pos[i-1];
205
206 // Expanding A.
207
208 if (AsB.uplo == 'U') {
209
210 for (i=n-1; i>=0; i--) {
211 pos[i] = index[i]+pcol[i+1]-pcol[i];
212 k = pos[i]-1;
213 for (j=pcol[i+1]-1; j>=pcol[i]; j--) {
214 value[k] = value[j];
215 irow[k--] = irow[j];
216 }
217 }
218 for (i=1; i<n; i++) {
219 k = index[i]+pcol[i+1]-pcol[i];
220 if ((k>index[i])&&(irow[k-1]==i)) k--;
221 for (j=index[i]; j<k; j++) {
222 value[pos[irow[j]]] = value[j];
223 irow[pos[irow[j]]++] = i;
224 }
225 }
226
227 }
228 else { // uplo == 'L'
229
230 for (i=n-1; i>=0; i--) {
231 k = index[i+1]-1;
232 for (j=pcol[i+1]-1; j>=pcol[i]; j--) {
233 value[k] = value[j];
234 irow[k--] = irow[j];
235 }
236 pos[i] = index[i];
237 }
238 for (i=0; i<(n-1); i++) {
239 k = index[i+1]-pcol[i+1]+pcol[i];
240 if ((k<index[i+1])&&(irow[k]==i)) k++;
241 for (j=k; j<index[i+1]; j++) {
242 value[pos[irow[j]]] = value[j];
243 irow[pos[irow[j]]++] = i;
244 }
245 }
246
247 }
248
249 AsB.nnz = index[n];
250
251 // Deleting temporary vectors.
252
253 delete[] pcol;
254 delete[] pos;
255
256 } // ExpandAsB.
257
258
259 template<class ARTYPE>
260 void ARchSymPencil<ARTYPE>::SubtractAsB(ARTYPE sigma)
261 {
262
263 int i, acol, bcol, asbcol, scol;
264
265 // Quitting function if A->uplo is not equal to B->uplo.
266
267 if ((A->uplo != B->uplo)&&(sigma != (ARTYPE)0)) {
268 throw ArpackError(ArpackError::DIFFERENT_TRIANGLES,
269 "ARchSymPencil::SubtractAsB");
270 }
271 AsB.uplo = A->uplo;
272
273 // Subtracting sigma*B from A.
274
275 AsB.index[0] = 0;
276 asbcol = 0;
277
278 for (i=0; i!=AsB.n; i++) {
279 bcol = B->pcol[i];
280 acol = A->pcol[i];
281 SparseSaxpy(-sigma, &B->a[bcol], &B->irow[bcol], B->pcol[i+1]-bcol,
282 &A->a[acol], &A->irow[acol], A->pcol[i+1]-acol,
283 &AsB.value[asbcol], &AsB.index[asbcol+AsB.n+1], scol);
284 asbcol += scol;
285 AsB.index[i+1] = asbcol;
286 }
287
288 // Expanding AsB.
289
290 ExpandAsB();
291
292 // Adding one to all elements of vector index
293 // because the decomposition function was written in FORTRAN.
294
295 for (i=0; i<=AsB.n+AsB.nnz; i++) AsB.index[i]++;
296
297 } // SubtractAsB.
298
299 */
300
301 template<class ARTYPE>
FactorAsB(ARTYPE sigma)302 void ARchSymPencil<ARTYPE>::FactorAsB(ARTYPE sigma)
303 {
304
305 // Quitting the function if A and B were not defined.
306
307 if (!(A->IsDefined()&&B->IsDefined())) {
308 throw ArpackError(ArpackError::DATA_UNDEFINED,
309 "ARchSymPencil::FactorAsB");
310 }
311
312
313 if (LAsB) cholmod_free_factor(&LAsB,&c);
314
315 cholmod_sparse* AsB;
316 if (sigma != 0.0)
317 {
318 std::cout << " Subtracting sigma B (sigma="<<sigma<<")"<<std::endl;
319 double alpha[2]; alpha[0]=1.0; alpha[1] = 1.0;
320 double beta[2]; beta[0] = -sigma; beta[1]=1.0;
321 AsB = cholmod_add(A->A,B->A,alpha,beta,1,0,&c);
322 }
323 else
324 AsB = A->A;
325
326 //FILE *fp;
327 //fp=fopen("AsB.asc", "w");
328 //cholmod_write_sparse(fp,AsB,NULL,NULL,&c);
329 //FILE *fpa;
330 //fpa=fopen("As.asc", "w");
331 //cholmod_write_sparse(fpa,B->A,NULL,NULL,&c);
332 //FILE *fpb;
333 //fpb=fopen("Bs.asc", "w");
334 //cholmod_write_sparse(fpb,A->A,NULL,NULL,&c);
335
336 LAsB = cholmod_analyze (AsB, &c) ;
337 int info = cholmod_factorize (AsB, LAsB, &c) ;
338
339 factoredAsB = (info != 0);
340 if (c.status != CHOLMOD_OK)
341 {
342 //std::cout << " sigma : " << sigma << std::endl;
343
344 Write_Cholmod_Sparse_Matrix("AsB-error.asc",AsB,&c);
345
346 throw ArpackError(ArpackError::INCONSISTENT_DATA,
347 "ARchSymPencil::FactorAsB");
348
349 factoredAsB = false;
350 }
351
352 if (sigma != 0.0)
353 cholmod_free_sparse(&AsB,&c);
354
355
356 } // FactorAsB (ARTYPE shift).
357
358
359 template<class ARTYPE>
MultInvBAv(ARTYPE * v,ARTYPE * w)360 void ARchSymPencil<ARTYPE>::MultInvBAv(ARTYPE* v, ARTYPE* w)
361 {
362
363 if (!B->IsFactored()) B->FactorA();
364
365 A->MultMv(v, w);
366 ::copy(A->ncols(), w, 1, v, 1);
367 B->MultInvv(w, w);
368
369 } // MultInvBAv.
370
371 template<class ARTYPE>
MultInvAsBv(ARTYPE * v,ARTYPE * w)372 void ARchSymPencil<ARTYPE>::MultInvAsBv(ARTYPE* v, ARTYPE* w)
373 {
374 if (!IsFactored()) {
375 throw ArpackError(ArpackError::NOT_FACTORED_MATRIX,
376 "ARchSymPencil::MultInvAsBv");
377 }
378
379 // Solving A.w = v (or AsI.w = v).
380
381 //create b from v (data is not copied!!)
382 cholmod_dense * b = Create_Cholmod_Dense_Matrix(A->n,1,v,&c);
383
384 cholmod_dense *x = cholmod_solve (CHOLMOD_A, LAsB, b, &c) ;
385
386 Get_Cholmod_Dense_Data(x, A->n, w);
387
388 free(b);
389 cholmod_free_dense(&x,&c);
390
391
392 } // MultInvAsBv
393
394 template<class ARTYPE>
395 inline void ARchSymPencil<ARTYPE>::
DefineMatrices(ARchSymMatrix<ARTYPE> & Ap,ARchSymMatrix<ARTYPE> & Bp)396 DefineMatrices(ARchSymMatrix<ARTYPE>& Ap, ARchSymMatrix<ARTYPE>& Bp)
397 {
398
399 A = &Ap;
400 B = &Bp;
401
402 if (A->n != B->n) {
403 throw ArpackError(ArpackError::INCOMPATIBLE_SIZES,
404 "ARchSymMatrix::DefineMatrices");
405 }
406
407 } // DefineMatrices.
408
409
410 template<class ARTYPE>
411 inline ARchSymPencil<ARTYPE>::
ARchSymPencil(ARchSymMatrix<ARTYPE> & Ap,ARchSymMatrix<ARTYPE> & Bp)412 ARchSymPencil(ARchSymMatrix<ARTYPE>& Ap, ARchSymMatrix<ARTYPE>& Bp)
413 {
414 cholmod_start (&c) ;
415 LAsB=NULL;
416 DefineMatrices(Ap, Bp);
417
418 } // Long constructor.
419
420
421 template<class ARTYPE>
422 ARchSymPencil<ARTYPE>& ARchSymPencil<ARTYPE>::
423 operator=(const ARchSymPencil<ARTYPE>& other)
424 {
425
426 if (this != &other) { // Stroustrup suggestion.
427 Copy(other);
428 }
429 return *this;
430
431 } // operator=.
432
433
434 #endif // ARUSPEN_H
435