1 //
2 // distsymm.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #include <iostream>
29 #include <math.h>
30
31 #include <util/misc/formio.h>
32 #include <util/keyval/keyval.h>
33 #include <math/scmat/dist.h>
34 #include <math/scmat/cmatrix.h>
35 #include <math/scmat/elemop.h>
36 #include <math/scmat/disthql.h>
37
38 using namespace std;
39 using namespace sc;
40
41 extern "C" { int DBmalloc_chain_check(const char *, int, int); }
42
43 /////////////////////////////////////////////////////////////////////////////
44 // DistSymmSCMatrix member functions
45
46 static ClassDesc DistSymmSCMatrix_cd(
47 typeid(DistSymmSCMatrix),"DistSymmSCMatrix",1,"public SymmSCMatrix",
48 0, 0, 0);
49
DistSymmSCMatrix(const RefSCDimension & a,DistSCMatrixKit * k)50 DistSymmSCMatrix::DistSymmSCMatrix(const RefSCDimension&a,DistSCMatrixKit*k):
51 SymmSCMatrix(a,k)
52 {
53 init_blocklist();
54 }
55
56 int
block_to_node(int i,int j) const57 DistSymmSCMatrix::block_to_node(int i, int j) const
58 {
59 if (j>i) {
60 ExEnv::errn() << indent << "DistSymmSCMatrix::block_to_node: j>i" << endl;
61 abort();
62 }
63
64 return ((i*(i+1))/2 + j)%messagegrp()->n();
65 }
66
67 Ref<SCMatrixBlock>
block_to_block(int i,int j) const68 DistSymmSCMatrix::block_to_block(int i, int j) const
69 {
70 if (j>i) {
71 ExEnv::errn() << indent << "DistSymmSCMatrix::block_to_block: j>i" << endl;
72 abort();
73 }
74
75 int offset = (i*(i+1))/2 + j;
76 int nproc = messagegrp()->n();
77
78 if ((offset%nproc) != messagegrp()->me()) return 0;
79
80 SCMatrixBlockListIter I;
81 for (I=blocklist->begin(); I!=blocklist->end(); I++) {
82 if (I.block()->blocki == i && I.block()->blockj == j)
83 return I.block();
84 }
85
86 ExEnv::errn() << indent << "DistSymmSCMatrix::block_to_block: internal error" << endl;
87 abort();
88 return 0;
89 }
90
91 double *
find_element(int i,int j) const92 DistSymmSCMatrix::find_element(int i, int j) const
93 {
94 if (j>i) {
95 int tmp = i; i=j; j=tmp;
96 }
97
98 int bi, oi;
99 d->blocks()->elem_to_block(i, bi, oi);
100
101 int bj, oj;
102 d->blocks()->elem_to_block(j, bj, oj);
103
104 Ref<SCMatrixBlock> ablk = block_to_block(bi, bj);
105 if (ablk.nonnull()) {
106 if (bi != bj) {
107 Ref<SCMatrixRectBlock> blk
108 = dynamic_cast<SCMatrixRectBlock*>(ablk.pointer());
109 if (blk.null()) return 0;
110 return &blk->data[oi*(blk->jend-blk->jstart)+oj];
111 }
112 else {
113 Ref<SCMatrixLTriBlock> blk
114 = dynamic_cast<SCMatrixLTriBlock*>(ablk.pointer());
115 if (blk.null()) return 0;
116 return &blk->data[(oi*(oi+1))/2+oj];
117 }
118 }
119 return 0;
120 }
121
122 int
element_to_node(int i,int j) const123 DistSymmSCMatrix::element_to_node(int i, int j) const
124 {
125 if (j>i) {
126 int tmp = i; i=j; j=tmp;
127 }
128
129 int bi, oi;
130 d->blocks()->elem_to_block(i, bi, oi);
131
132 int bj, oj;
133 d->blocks()->elem_to_block(j, bj, oj);
134
135 return block_to_node(bi,bj);
136 }
137
138 void
init_blocklist()139 DistSymmSCMatrix::init_blocklist()
140 {
141 int i, j, index;
142 int nproc = messagegrp()->n();
143 int me = messagegrp()->me();
144 SCMatrixBlock *b;
145 blocklist = new SCMatrixBlockList;
146 for (i=0, index=0; i<d->blocks()->nblock(); i++) {
147 for (j=0; j<i; j++, index++) {
148 if (index%nproc != me) continue;
149 b = new SCMatrixRectBlock(d->blocks()->start(i),
150 d->blocks()->fence(i),
151 d->blocks()->start(j),
152 d->blocks()->fence(j));
153 b->blocki = i;
154 b->blockj = j;
155 blocklist->insert(b);
156 }
157 if (index%nproc == me) {
158 b = new SCMatrixLTriBlock(d->blocks()->start(i),
159 d->blocks()->fence(i));
160 b->blocki = i;
161 b->blockj = i;
162 blocklist->insert(b);
163 }
164 index++;
165 }
166 }
167
~DistSymmSCMatrix()168 DistSymmSCMatrix::~DistSymmSCMatrix()
169 {
170 }
171
172 double
get_element(int i,int j) const173 DistSymmSCMatrix::get_element(int i,int j) const
174 {
175 double res;
176 double *e = find_element(i,j);
177 if (e) {
178 res = *e;
179 messagegrp()->bcast(res, messagegrp()->me());
180 }
181 else {
182 messagegrp()->bcast(res, element_to_node(i, j));
183 }
184 return res;
185 }
186
187 void
set_element(int i,int j,double a)188 DistSymmSCMatrix::set_element(int i,int j,double a)
189 {
190 double *e = find_element(i,j);
191 if (e) {
192 *e = a;
193 }
194 }
195
196 void
accumulate_element(int i,int j,double a)197 DistSymmSCMatrix::accumulate_element(int i,int j,double a)
198 {
199 double *e = find_element(i,j);
200 if (e) {
201 *e += a;
202 }
203 }
204
205 SymmSCMatrix *
get_subblock(int br,int er)206 DistSymmSCMatrix::get_subblock(int br, int er)
207 {
208 error("get_subblock");
209 return 0;
210 }
211
212 SCMatrix *
get_subblock(int,int,int,int)213 DistSymmSCMatrix::get_subblock(int, int, int, int)
214 {
215 error("get_subblock");
216 return 0;
217 }
218
219 void
assign_subblock(SCMatrix * sb,int br,int er,int bc,int ec)220 DistSymmSCMatrix::assign_subblock(SCMatrix*sb, int br, int er, int bc, int ec)
221 {
222 error("assign_subblock");
223 }
224
225 void
assign_subblock(SymmSCMatrix * sb,int br,int er)226 DistSymmSCMatrix::assign_subblock(SymmSCMatrix*sb, int br, int er)
227 {
228 error("accumulate_subblock");
229 }
230
231 void
accumulate_subblock(SCMatrix * sb,int br,int er,int bc,int ec)232 DistSymmSCMatrix::accumulate_subblock(SCMatrix*sb, int br, int er, int bc, int ec)
233 {
234 error("accumulate_subblock");
235 }
236
237 void
accumulate_subblock(SymmSCMatrix * sb,int br,int er)238 DistSymmSCMatrix::accumulate_subblock(SymmSCMatrix*sb, int br, int er)
239 {
240 error("accumulate_subblock");
241 }
242
243 SCVector *
get_row(int i)244 DistSymmSCMatrix::get_row(int i)
245 {
246 error("get_row");
247 return 0;
248 }
249
250 void
assign_row(SCVector * v,int i)251 DistSymmSCMatrix::assign_row(SCVector *v, int i)
252 {
253 error("assign_row");
254 }
255
256 void
accumulate_row(SCVector * v,int i)257 DistSymmSCMatrix::accumulate_row(SCVector *v, int i)
258 {
259 error("accumulate_row");
260 }
261
262 void
accumulate(const SymmSCMatrix * a)263 DistSymmSCMatrix::accumulate(const SymmSCMatrix*a)
264 {
265 // make sure that the arguments is of the correct type
266 const DistSymmSCMatrix* la
267 = require_dynamic_cast<const DistSymmSCMatrix*>(a,"DistSymmSCMatrix::accumulate");
268
269 // make sure that the dimensions match
270 if (!dim()->equiv(la->dim())) {
271 ExEnv::errn() << indent << "DistSymmSCMatrix::accumulate(SCMatrix*a): "
272 << "dimensions don't match\n";
273 abort();
274 }
275
276 SCMatrixBlockListIter i1, i2;
277 for (i1=la->blocklist->begin(),i2=blocklist->begin();
278 i1!=la->blocklist->end() && i2!=blocklist->end();
279 i1++,i2++) {
280 int n = i1.block()->ndat();
281 if (n != i2.block()->ndat()) {
282 ExEnv::errn() << indent
283 << "DistSymmSCMatrixListSubblockIter::accumulate block "
284 << "mismatch: internal error" << endl;
285 abort();
286 }
287 double *dat1 = i1.block()->dat();
288 double *dat2 = i2.block()->dat();
289 for (int i=0; i<n; i++) {
290 dat2[i] += dat1[i];
291 }
292 }
293 }
294
295 double
invert_this()296 DistSymmSCMatrix::invert_this()
297 {
298 RefDiagSCMatrix refa = kit()->diagmatrix(d);
299 RefSCMatrix refb = kit()->matrix(d,d);
300 diagonalize(refa.pointer(),refb.pointer());
301 double determ = 1.0;
302 for (int i=0; i<dim()->n(); i++) {
303 double val = refa->get_element(i);
304 determ *= val;
305 }
306 Ref<SCElementOp> op = new SCElementInvert(1.0e-12);
307 refa->element_op(op.pointer());
308 assign(0.0);
309 accumulate_transform(refb.pointer(), refa.pointer());
310 return determ;
311 }
312
313 double
determ_this()314 DistSymmSCMatrix::determ_this()
315 {
316 return invert_this();
317 }
318
319 double
trace()320 DistSymmSCMatrix::trace()
321 {
322 double ret=0.0;
323 Ref<SCMatrixSubblockIter> I = local_blocks(SCMatrixSubblockIter::Read);
324 for (I->begin(); I->ready(); I->next()) {
325 Ref<SCMatrixLTriBlock> b = dynamic_cast<SCMatrixLTriBlock*>(I->block());
326 if (b.nonnull() && b->blocki == b->blockj) {
327 int ni = b->end-b->start;
328 double *data = b->data;
329 for (int i=0; i<ni; i++) {
330 data += i;
331 ret += *data;
332 }
333 }
334 }
335 messagegrp()->sum(ret);
336
337 return ret;
338 }
339
340 double
solve_this(SCVector * v)341 DistSymmSCMatrix::solve_this(SCVector*v)
342 {
343 DistSCVector* lv =
344 require_dynamic_cast<DistSCVector*>(v,"DistSymmSCMatrix::solve_this");
345
346 // make sure that the dimensions match
347 if (!dim()->equiv(lv->dim())) {
348 ExEnv::errn() << indent << "DistSymmSCMatrix::solve_this(SCVector*v): "
349 << "dimensions don't match\n";
350 abort();
351 }
352
353 error("no solve this");
354
355 return 0.0;
356 }
357
358 void
gen_invert_this()359 DistSymmSCMatrix::gen_invert_this()
360 {
361 invert_this();
362 }
363
364 void
diagonalize(DiagSCMatrix * a,SCMatrix * b)365 DistSymmSCMatrix::diagonalize(DiagSCMatrix*a,SCMatrix*b)
366 {
367 const char* name = "DistSymmSCMatrix::diagonalize";
368 // make sure that the argument are of the correct type
369 require_dynamic_cast<DistDiagSCMatrix*>(a,name);
370 DistSCMatrix* lb = require_dynamic_cast<DistSCMatrix*>(b,name);
371
372 int n = dim()->n();
373 int me = messagegrp()->me();
374 int nproc = messagegrp()->n();
375
376 RefSCMatrix arect = kit()->matrix(dim(),dim());
377 DistSCMatrix *rect = dynamic_cast<DistSCMatrix*>(arect.pointer());
378 rect->assign(0.0);
379 rect->accumulate(this);
380
381 // This sets up the index list of columns to be stored on this node
382 int nvec = n/nproc + (me<(n%nproc)?1:0);
383 int *ivec = new int[nvec];
384 for (int i=0; i<nvec; i++) {
385 ivec[i] = i*nproc + me;
386 }
387
388 rect->create_vecform(DistSCMatrix::Col,nvec);
389 rect->vecform_op(DistSCMatrix::CopyToVec,ivec);
390 lb->create_vecform(DistSCMatrix::Col,nvec);
391
392 double *d = new double[n];
393 dist_diagonalize(n, rect->nvec, rect->vec[0], d, lb->vec[0],
394 messagegrp());
395
396 // put d into the diagonal matrix
397 a->assign(d);
398
399 lb->vecform_op(DistSCMatrix::CopyFromVec, ivec);
400 lb->delete_vecform();
401 rect->delete_vecform();
402 arect = 0;
403 delete[] ivec;
404 }
405
406 // computes this += a + a.t
407 void
accumulate_symmetric_sum(SCMatrix * a)408 DistSymmSCMatrix::accumulate_symmetric_sum(SCMatrix*a)
409 {
410 // make sure that the argument is of the correct type
411 DistSCMatrix* la
412 = require_dynamic_cast<DistSCMatrix*>(a,"DistSymmSCMatrix::"
413 "accumulate_symmetric_sum");
414
415 if (!dim()->equiv(la->rowdim()) || !dim()->equiv(la->coldim())) {
416 ExEnv::errn() << indent << "DistSymmSCMatrix::"
417 << "accumulate_symmetric_sum(SCMatrix*a): bad dim\n";
418 abort();
419 }
420
421 Ref<SCMatrixSubblockIter> I = all_blocks(SCMatrixSubblockIter::Accum);
422 for (I->begin(); I->ready(); I->next()) {
423 Ref<SCMatrixBlock> block = I->block();
424 // see if i've got this block
425 Ref<SCMatrixBlock> localblock
426 = la->block_to_block(block->blocki,block->blockj);
427 if (localblock.nonnull()) {
428 // the diagonal blocks require special handling
429 if (block->blocki == block->blockj) {
430 int n = la->rowblocks()->size(block->blocki);
431 double *dat1 = block->dat();
432 double *dat2 = localblock->dat();
433 for (int i=0; i<n; i++) {
434 for (int j=0; j<=i; j++) {
435 double tmp = 0.0;
436 tmp += dat2[i*n+j];
437 tmp += dat2[j*n+i];
438 *dat1 += tmp;
439 dat1++;
440 }
441 }
442 }
443 else {
444 int n = block->ndat();
445 double *dat1 = block->dat();
446 double *dat2 = localblock->dat();
447 for (int i=0; i<n; i++) {
448 dat1[i] += dat2[i];
449 }
450 }
451 }
452 // now for the transpose
453 if (block->blocki != block->blockj) {
454 localblock = la->block_to_block(block->blockj,block->blocki);
455 if (localblock.nonnull()) {
456 int nr = la->rowblocks()->size(block->blocki);
457 int nc = la->rowblocks()->size(block->blockj);
458 double *dat1 = block->dat();
459 double *dat2 = localblock->dat();
460 for (int i=0; i<nr; i++) {
461 for (int j=0; j<nc; j++) {
462 *dat1++ += dat2[j*nr+i];
463 }
464 }
465 }
466 }
467 }
468 }
469
470 void
element_op(const Ref<SCElementOp> & op)471 DistSymmSCMatrix::element_op(const Ref<SCElementOp>& op)
472 {
473 SCMatrixBlockListIter i;
474 for (i = blocklist->begin(); i != blocklist->end(); i++) {
475 op->process_base(i.block());
476 }
477 if (op->has_collect()) op->collect(messagegrp());
478 }
479
480 void
element_op(const Ref<SCElementOp2> & op,SymmSCMatrix * m)481 DistSymmSCMatrix::element_op(const Ref<SCElementOp2>& op,
482 SymmSCMatrix* m)
483 {
484 DistSymmSCMatrix *lm
485 = require_dynamic_cast<DistSymmSCMatrix*>(m,"DistSymSCMatrix::element_op");
486
487 if (!dim()->equiv(lm->dim())) {
488 ExEnv::errn() << indent << "DistSymmSCMatrix: bad element_op\n";
489 abort();
490 }
491 SCMatrixBlockListIter i, j;
492 for (i = blocklist->begin(), j = lm->blocklist->begin();
493 i != blocklist->end();
494 i++, j++) {
495 op->process_base(i.block(), j.block());
496 }
497 if (op->has_collect()) op->collect(messagegrp());
498 }
499
500 void
element_op(const Ref<SCElementOp3> & op,SymmSCMatrix * m,SymmSCMatrix * n)501 DistSymmSCMatrix::element_op(const Ref<SCElementOp3>& op,
502 SymmSCMatrix* m,SymmSCMatrix* n)
503 {
504 DistSymmSCMatrix *lm
505 = require_dynamic_cast<DistSymmSCMatrix*>(m,"DistSymSCMatrix::element_op");
506 DistSymmSCMatrix *ln
507 = require_dynamic_cast<DistSymmSCMatrix*>(n,"DistSymSCMatrix::element_op");
508
509 if (!dim()->equiv(lm->dim()) || !dim()->equiv(ln->dim())) {
510 ExEnv::errn() << indent << "DistSymmSCMatrix: bad element_op\n";
511 abort();
512 }
513 SCMatrixBlockListIter i, j, k;
514 for (i = blocklist->begin(),
515 j = lm->blocklist->begin(),
516 k = ln->blocklist->begin();
517 i != blocklist->end();
518 i++, j++, k++) {
519 op->process_base(i.block(), j.block(), k.block());
520 }
521 if (op->has_collect()) op->collect(messagegrp());
522 }
523
524 Ref<SCMatrixSubblockIter>
local_blocks(SCMatrixSubblockIter::Access access)525 DistSymmSCMatrix::local_blocks(SCMatrixSubblockIter::Access access)
526 {
527 return new SCMatrixListSubblockIter(access, blocklist);
528 }
529
530 Ref<SCMatrixSubblockIter>
all_blocks(SCMatrixSubblockIter::Access access)531 DistSymmSCMatrix::all_blocks(SCMatrixSubblockIter::Access access)
532 {
533 return new DistSCMatrixListSubblockIter(access, blocklist, messagegrp());
534 }
535
536 void
error(const char * msg)537 DistSymmSCMatrix::error(const char *msg)
538 {
539 ExEnv::errn() << "DistSymmSCMatrix: error: " << msg << endl;
540 }
541
542 Ref<DistSCMatrixKit>
skit()543 DistSymmSCMatrix::skit()
544 {
545 return dynamic_cast<DistSCMatrixKit*>(kit().pointer());
546 }
547
548 void
convert_accumulate(SymmSCMatrix * a)549 DistSymmSCMatrix::convert_accumulate(SymmSCMatrix*a)
550 {
551 SymmSCMatrix::convert_accumulate(a);
552
553 #if 0
554 DistSymmSCMatrix *d = require_dynamic_cast<DistSymmSCMatrix*>(a,
555 "DistSymmSCMatrix::convert_accumulate");
556
557 SCMatrixBlockListIter i, j;
558 for (i = blocklist->begin(), j = d->blocklist->begin();
559 i != blocklist->end();
560 i++, j++) {
561 }
562 #endif
563 }
564
565 /////////////////////////////////////////////////////////////////////////////
566
567 // Local Variables:
568 // mode: c++
569 // c-file-style: "CLJ"
570 // End:
571