/[ascend]/trunk/ascend4/solver/mtx_query.c
ViewVC logotype

Contents of /trunk/ascend4/solver/mtx_query.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 50034 byte(s)
Setting up web subdirectory in repository
1 /*
2 * mtx: Ascend Sparse Matrix Package
3 * by Karl Michael Westerberg
4 * Created: 2/6/90
5 * Version: $Revision: 1.13 $
6 * Version control file: $RCSfile: mtx_query.c,v $
7 * Date last modified: $Date: 1997/07/28 20:53:07 $
8 * Last modified by: $Author: ballan $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1990 Karl Michael Westerberg
13 * Copyright (C) 1993 Joseph Zaher
14 * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15 * Copyright (C) 1995 Kirk Andre Abbott, Benjamin Andrew Allan
16 * Copyright (C) 1996 Benjamin Andrew Allan
17 *
18 * The SLV solver is free software; you can redistribute
19 * it and/or modify it under the terms of the GNU General Public License as
20 * published by the Free Software Foundation; either version 2 of the
21 * License, or (at your option) any later version.
22 *
23 * The SLV solver is distributed in hope that it will be
24 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
25 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
26 * General Public License for more details.
27 *
28 * You should have received a copy of the GNU General Public License along with
29 * the program; if not, write to the Free Software Foundation, Inc., 675
30 * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
31 * COPYING is found in ../compiler.
32 */
33
34 #include <math.h>
35 #include "utilities/ascConfig.h"
36 #include "utilities/ascMalloc.h"
37 #include "utilities/mem.h"
38 #include "solver/mtx.h"
39 /* grab our private parts */
40 #define __MTX_C_SEEN__
41 #include "solver/mtx_use_only.h"
42
43 real64 mtx_next_in_row( mtx_matrix_t mtx, mtx_coord_t *coord, mtx_range_t *rng)
44 {
45 static struct element_t *elt = NULL;
46 struct element_t Rewind;
47
48 #if MTX_DEBUG
49 if(!mtx_check_matrix(mtx)) return D_ZERO;
50 #endif
51 if( coord->col == mtx_FIRST ) {
52 if (not_in_range(0,mtx->order-1,coord->row)) {
53 /* coord->col = mtx_LAST; already true */
54 return D_ZERO;
55 }
56 Rewind.next.col = mtx->hdr.row[mtx->perm.row.cur_to_org[coord->row]];
57 elt = mtx_next_col(&Rewind,rng,mtx->perm.col.org_to_cur);
58 } else elt = mtx_next_col(elt,rng,mtx->perm.col.org_to_cur);
59 coord->col = (ISNULL(elt)) ? mtx_LAST : mtx->perm.col.org_to_cur[elt->col];
60 return( ISNULL(elt) ? D_ZERO : elt->value );
61 }
62
63 real64 mtx_next_in_col( mtx_matrix_t mtx, mtx_coord_t *coord, mtx_range_t *rng)
64 {
65 static struct element_t *elt = NULL;
66 struct element_t Rewind;
67
68 #if MTX_DEBUG
69 if(!mtx_check_matrix(mtx)) return 0;
70 #endif
71 if( coord->row == mtx_FIRST ) {
72 if (not_in_range(0,mtx->order-1,coord->col)) {
73 /* coord->row = mtx_LAST; already true */
74 return D_ZERO;
75 }
76 Rewind.next.row = mtx->hdr.col[mtx->perm.col.cur_to_org[coord->col]];
77 elt = mtx_next_row(&Rewind,rng,mtx->perm.row.org_to_cur);
78 } else elt = mtx_next_row(elt,rng,mtx->perm.row.org_to_cur);
79 coord->row = (ISNULL(elt)) ? mtx_LAST : mtx->perm.row.org_to_cur[elt->row];
80 return( ISNULL(elt) ? D_ZERO : elt->value );
81 }
82
83 real64 mtx_row_max( mtx_matrix_t mtx, mtx_coord_t *coord, mtx_range_t *rng, real64 *val)
84 {
85 struct element_t *elt,*maxelt;
86 real64 maxabs;
87 int32 *tocur;
88
89 #if MTX_DEBUG
90 if(!mtx_check_matrix(mtx)) return (double)-1.0;
91 #endif
92 maxabs = D_ZERO;
93 maxelt = NULL;
94 tocur = mtx->perm.col.org_to_cur;
95 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[coord->row]];
96 if( rng == mtx_ALL_COLS ) {
97 for( ; NOTNULL(elt); elt = elt->next.col ) {
98 if( fabs(elt->value) < maxabs ) {
99 continue; /* skip uninteresting elements */
100 }
101 if (ISNULL(maxelt)) { /* first element found */
102 maxelt = elt;
103 maxabs = fabs(elt->value);
104 } else { /* tie breakers: elt >= maxelt */
105 if( fabs(elt->value) > maxabs /* bigger */ ||
106 tocur[elt->col] < tocur[maxelt->col] /* nearer */ ) {
107 maxelt = elt;
108 maxabs = fabs(elt->value);
109 }
110 }
111 }
112 } else {
113 if( rng->high >= rng->low ) {
114 for( ; NOTNULL(elt); elt = elt->next.col ) {
115 if( in_range(rng,tocur[elt->col]) ) {
116 if( fabs(elt->value) < maxabs ) {
117 continue; /* skip uninteresting elements */
118 }
119 if (ISNULL(maxelt)) { /* first element found */
120 maxelt = elt;
121 maxabs = fabs(elt->value);
122 } else { /* tie breakers: elt >= maxelt */
123 if( fabs(elt->value) > maxabs /* bigger */ ||
124 tocur[elt->col] < tocur[maxelt->col] /* nearer */ ) {
125 maxelt = elt;
126 maxabs = fabs(elt->value);
127 }
128 }
129 }
130 }
131 }
132 }
133
134 coord->col = (ISNULL(maxelt)) ? mtx_NONE : tocur[maxelt->col];
135 if (val) {
136 if (maxelt) *val=maxelt->value; else *val=D_ZERO;
137 }
138 return(maxabs);
139 }
140
141 real64 mtx_col_max( mtx_matrix_t mtx, mtx_coord_t *coord, mtx_range_t *rng, real64 *val)
142 {
143 struct element_t *elt,*maxelt;
144 real64 maxabs;
145 int32 *tocur;
146
147 #if MTX_DEBUG
148 if(!mtx_check_matrix(mtx)) return (double)-1.0;
149 #endif
150 maxabs = D_ZERO;
151 maxelt = NULL;
152 tocur = mtx->perm.row.org_to_cur;
153 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[coord->col]];
154 if( rng == mtx_ALL_ROWS ) {
155 for( ; NOTNULL(elt); elt = elt->next.row ) {
156 if( fabs(elt->value) < maxabs ) {
157 continue; /* skip uninteresting elements */
158 }
159 if (ISNULL(maxelt)) { /* first element found */
160 maxelt = elt;
161 maxabs = fabs(elt->value);
162 } else { /* tie breakers: elt >= maxelt */
163 if( fabs(elt->value) > maxabs /* bigger */ ||
164 tocur[elt->row] < tocur[maxelt->row] /* nearer */ ) {
165 maxelt = elt;
166 maxabs = fabs(elt->value);
167 }
168 }
169 }
170 } else {
171 if( rng->high >= rng->low ) {
172 for( ; NOTNULL(elt); elt = elt->next.row ) {
173 if( in_range(rng,tocur[elt->row]) ) {
174 if( fabs(elt->value) < maxabs ) {
175 continue; /* skip uninteresting elements */
176 }
177 if (ISNULL(maxelt)) { /* first element found */
178 maxelt = elt;
179 maxabs = fabs(elt->value);
180 } else { /* tie breakers: elt >= maxelt */
181 if( fabs(elt->value) > maxabs /* bigger */ ||
182 tocur[elt->row] < tocur[maxelt->row] /* nearer */ ) {
183 maxelt = elt;
184 maxabs = fabs(elt->value);
185 }
186 }
187 }
188 }
189 }
190 }
191
192 coord->row = (ISNULL(maxelt)) ? mtx_NONE : tocur[maxelt->row];
193 if (val) {
194 if (maxelt) *val=maxelt->value; else *val=D_ZERO;
195 }
196 return(maxabs);
197 }
198
199 real64 mtx_row_min( mtx_matrix_t mtx, mtx_coord_t *coord,
200 mtx_range_t *rng, real64 *val, real64 minval)
201 {
202 struct element_t *elt,*minelt;
203 real64 minabs;
204 real64 big_num = 1e50;
205 int32 *tocur;
206
207 #if MTX_DEBUG
208 if(!mtx_check_matrix(mtx)) return (double)-1.0;
209 #endif
210 minabs = big_num;
211 minelt = NULL;
212 tocur = mtx->perm.col.org_to_cur;
213 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[coord->row]];
214 if( rng == mtx_ALL_COLS ) {
215 for( ; NOTNULL(elt); elt = elt->next.col ) {
216 if( fabs(elt->value) > minabs ) {
217 continue; /* skip uninteresting elements */
218 }
219 if( (fabs(elt->value) < minabs /* smaller */ ||
220 tocur[elt->col] < tocur[minelt->col] /* nearer */ )
221 && (fabs(elt->value) > minval)) {
222 minelt = elt;
223 minabs = fabs(elt->value);
224 }
225 }
226 } else {
227 if( rng->high >= rng->low ) {
228 for( ; NOTNULL(elt); elt = elt->next.col ) {
229 if( in_range(rng,tocur[elt->col]) ) {
230 if( fabs(elt->value) > minabs ) {
231 continue; /* skip uninteresting elements */
232 }
233 if( (fabs(elt->value) < minabs /* smaller */ ||
234 tocur[elt->col] < tocur[minelt->col] /* nearer */ )
235 && (fabs(elt->value) > minval)) {
236 minelt = elt;
237 minabs = fabs(elt->value);
238 }
239 }
240 }
241 }
242 }
243 if(minabs == big_num){
244 coord->col = mtx_NONE;
245 if (val) {
246 *val=0;
247 }
248 minabs = 1;
249 } else {
250 coord->col = (ISNULL(minelt)) ? mtx_NONE : tocur[minelt->col];
251 if (val) {
252 *val=minelt->value;
253 }
254 }
255 return(minabs);
256 }
257
258 real64 mtx_col_min( mtx_matrix_t mtx, mtx_coord_t *coord,
259 mtx_range_t *rng, real64 *val, real64 minval)
260 {
261 struct element_t *elt,*minelt;
262 real64 minabs;
263 real64 big_num = 1e50;
264 int32 *tocur;
265
266 #if MTX_DEBUG
267 if(!mtx_check_matrix(mtx)) return (double)-1.0;
268 #endif
269 minabs = big_num;
270 minelt = NULL;
271 tocur = mtx->perm.row.org_to_cur;
272 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[coord->col]];
273 if( rng == mtx_ALL_ROWS ) {
274 for( ; NOTNULL(elt); elt = elt->next.row ) {
275 if( fabs(elt->value) > minabs ) {
276 continue; /* skip uninteresting elements */
277 }
278 if( (fabs(elt->value) < minabs /* smaller */ ||
279 tocur[elt->row] < tocur[minelt->row] /* nearer */ )
280 && (fabs(elt->value) > minval)) {
281 minelt = elt;
282 minabs = fabs(elt->value);
283 }
284 }
285 } else {
286 if( rng->high >= rng->low ) {
287 for( ; NOTNULL(elt); elt = elt->next.row ) {
288 if( in_range(rng,tocur[elt->row]) ) {
289 if( fabs(elt->value) > minabs ) {
290 continue; /* skip uninteresting elements */
291 }
292 if( (fabs(elt->value) < minabs /* smaller */ ||
293 tocur[elt->row] < tocur[minelt->row] /* nearer */ )
294 && (fabs(elt->value) > minval)) {
295 minelt = elt;
296 minabs = fabs(elt->value);
297 }
298 }
299 }
300 }
301 }
302 if(minabs == big_num){
303 coord->col = mtx_NONE;
304 if (val) {
305 *val=0;
306 }
307 minabs = 1;
308 } else {
309 coord->row = (ISNULL(minelt)) ? mtx_NONE : tocur[minelt->row];
310 if (val) {
311 *val=minelt->value;
312 }
313 }
314 return(minabs);
315 }
316
317 real64 mtx_get_pivot_col( mtx_matrix_t mtx, mtx_coord_t *coord,
318 mtx_range_t *rng, real64 *val,
319 real64 ptol, real64 eps)
320 {
321 struct element_t *elt;
322 real64 rmax, thold;
323 real64 *rvec;
324 int32 *tocur;
325 int32 *ivec;
326 char *mark;
327 int32 space,j, k,len,pivot;
328
329
330 #if MTX_DEBUG
331 if(!mtx_check_matrix(mtx)) return (double)-1.0;
332 #endif
333 if (rng == mtx_ALL_COLS) {
334 space = mtx->order;
335 } else {
336 space = rng->high - rng->low +1;
337 }
338 if (space <1 || coord->row < 0 || coord->row >= mtx->order) {
339 /* empty or inverted range is bogus */
340 coord->col = mtx_NONE;
341 return D_ZERO;
342 }
343 mark = mtx_null_mark(space);
344 if (ISNULL(mark)) {
345 FPRINTF(g_mtxerr,"Error: (mtx.c) mtx_get_pivot_col. Insufficient memory.\n");
346 coord->col = mtx_NONE;
347 return D_ZERO;
348 }
349 rvec = mtx_null_sum(space);
350 if (ISNULL(rvec)) {
351 FPRINTF(g_mtxerr,"Error: (mtx.c) mtx_get_pivot_col. Insufficient memory.\n");
352 mtx_null_mark_release();
353 coord->col = mtx_NONE;
354 return D_ZERO;
355 }
356 ivec = mtx_null_index(space);
357 if (ISNULL(ivec)) {
358 FPRINTF(g_mtxerr,"Error: (mtx.c) mtx_get_pivot_col. Insufficient memory.\n");
359 mtx_null_mark_release();
360 mtx_null_sum_release();
361 coord->col = mtx_NONE;
362 return D_ZERO;
363 }
364 /* we're happy now, init search */
365
366 rmax = D_ZERO;
367 thold = eps = fabs(eps);
368 len = 0;
369 tocur = mtx->perm.col.org_to_cur;
370 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[coord->row]];
371
372 if( rng == mtx_ALL_COLS ) {
373 for( ; NOTNULL(elt); elt = elt->next.col ) {
374 if( (rvec[len] = fabs(elt->value) ) > thold ) {
375 mark[len] = (elt->value < D_ZERO);
376 ivec[len] = tocur[elt->col];
377 if (rvec[len] > rmax) {
378 rmax = rvec[len];
379 thold = MAX(rvec[len]*ptol,eps);
380 }
381 len++;
382 }
383 }
384 if (!len) {
385 /* found nothing > eps */
386 coord->col = mtx_NONE;
387 rvec[len] = D_ZERO; /* rezero the only buff we've touched. */
388 mtx_null_mark_release();
389 mtx_null_sum_release();
390 mtx_null_index_release();
391 return D_ZERO;
392 }
393 k = mtx->order;
394 } else {
395 for( ; NOTNULL(elt); elt = elt->next.col ) {
396 if( in_range(rng,tocur[elt->col]) ) {
397 if( (rvec[len] = fabs(elt->value) ) > thold ) {
398 mark[len] = (elt->value < D_ZERO);
399 ivec[len] = tocur[elt->col];
400 if (rvec[len] > rmax) {
401 rmax = rvec[len];
402 thold = MAX(rvec[len]*ptol,eps);
403 }
404 len++;
405 }
406 }
407 }
408 if (!len) {
409 /* found nothing > eps */
410 coord->col = mtx_NONE;
411 rvec[len] = D_ZERO; /* rezero the only buff we've touched. */
412 mtx_null_mark_release();
413 mtx_null_sum_release();
414 mtx_null_index_release();
415 return D_ZERO;
416 }
417 k = rng->high+1;
418 }
419
420 /* thold is now ptol * maxabsinrow, and k 1 past the right edge of rng. */
421 /* k becomes the nz.col we're after, pivot is the rvec/ivec location of k */
422 for (j = 0, pivot = k; j < len; j++) {
423 if (rvec[j] >= thold && ivec[j] < k ) {
424 k = ivec[j];
425 pivot = j;
426 }
427 }
428 coord->col = k; /* set col return */
429 rmax = rvec[pivot]; /* set return number */
430 if (val) { /* set val return */
431 if (mark[pivot]) {
432 *val = -rvec[pivot];
433 } else {
434 *val = rvec[pivot];
435 }
436 }
437 /* clean up buffers */
438 mtx_zero_real64(rvec,len);
439 mtx_zero_int32(ivec,len);
440 mtx_zero_char(mark,len);
441 mtx_null_mark_release();
442 mtx_null_sum_release();
443 mtx_null_index_release();
444 /* quit */
445 return(rmax);
446 }
447
448 real64 mtx_get_pivot_row( mtx_matrix_t mtx, mtx_coord_t *coord,
449 mtx_range_t *rng, real64 *val,
450 real64 ptol, real64 eps)
451 {
452 struct element_t *elt;
453 real64 rmax, thold;
454 real64 *rvec;
455 int32 *tocur;
456 int32 *ivec;
457 char *mark;
458 int32 space,j, k,len,pivot;
459
460
461 #if MTX_DEBUG
462 if(!mtx_check_matrix(mtx)) return (double)-1.0;
463 #endif
464 if (rng == mtx_ALL_ROWS) {
465 space = mtx->order;
466 } else {
467 space = rng->high - rng->low +1;
468 }
469 if (space <1 || coord->col < 0 || coord->col >= mtx->order) {
470 /* empty or inverted range is bogus */
471 coord->row = mtx_NONE;
472 return D_ZERO;
473 }
474 mark = mtx_null_mark(space);
475 if (ISNULL(mark)) {
476 FPRINTF(g_mtxerr,"Error: (mtx.c) mtx_get_pivot_row. Insufficient memory.\n");
477 coord->row = mtx_NONE;
478 return D_ZERO;
479 }
480 rvec = mtx_null_sum(space);
481 if (ISNULL(rvec)) {
482 FPRINTF(g_mtxerr,"Error: (mtx.c) mtx_get_pivot_row. Insufficient memory.\n");
483 mtx_null_mark_release();
484 coord->row = mtx_NONE;
485 return D_ZERO;
486 }
487 ivec = mtx_null_index(space);
488 if (ISNULL(ivec)) {
489 FPRINTF(g_mtxerr,"Error: (mtx.c) mtx_get_pivot_row. Insufficient memory.\n");
490 mtx_null_mark_release();
491 mtx_null_sum_release();
492 coord->row = mtx_NONE;
493 return D_ZERO;
494 }
495 /* we're happy now, init search */
496
497 rmax = D_ZERO;
498 thold = eps = fabs(eps);
499 len = 0;
500 tocur = mtx->perm.row.org_to_cur;
501 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[coord->col]];
502
503 if( rng == mtx_ALL_ROWS ) {
504 for( ; NOTNULL(elt); elt = elt->next.row ) {
505 if( (rvec[len] = fabs(elt->value) ) > thold ) {
506 mark[len] = (elt->value < D_ZERO);
507 ivec[len] = tocur[elt->row];
508 if (rvec[len] > rmax) {
509 rmax = rvec[len];
510 thold = MAX(rvec[len]*ptol,eps);
511 }
512 len++;
513 }
514 }
515 if (!len) {
516 /* found nothing > eps */
517 coord->row = mtx_NONE;
518 rvec[len] = D_ZERO; /* rezero the only buff we've touched. */
519 mtx_null_mark_release();
520 mtx_null_sum_release();
521 mtx_null_index_release();
522 return D_ZERO;
523 }
524 k = mtx->order;
525 } else {
526 for( ; NOTNULL(elt); elt = elt->next.row ) {
527 if( in_range(rng,tocur[elt->row]) ) {
528 if( (rvec[len] = fabs(elt->value) ) > thold ) {
529 mark[len] = (elt->value < D_ZERO);
530 ivec[len] = tocur[elt->row];
531 if (rvec[len] > rmax) {
532 rmax = rvec[len];
533 thold = MAX(rvec[len]*ptol,eps);
534 }
535 len++;
536 }
537 }
538 }
539 if (!len) {
540 /* found nothing > eps */
541 coord->row = mtx_NONE;
542 rvec[len] = D_ZERO; /* rezero the only buff we've touched. */
543 mtx_null_mark_release();
544 mtx_null_sum_release();
545 mtx_null_index_release();
546 return D_ZERO;
547 }
548 k = rng->high+1;
549 }
550
551 /* thold is now ptol * maxabsincol, and k 1 past the bottom edge of rng. */
552 /* k becomes the nz.col we're after, pivot is the rvec/ivec location of k */
553 for (j = 0, pivot = k; j < len; j++) {
554 if (rvec[j] >= thold && ivec[j] < k ) {
555 k = ivec[j];
556 pivot = j;
557 }
558 }
559 coord->row = k; /* set row return */
560 rmax = rvec[pivot]; /* set return number */
561 if (val) { /* set val return */
562 if (mark[pivot]) {
563 *val = -rvec[pivot];
564 } else {
565 *val = rvec[pivot];
566 }
567 }
568 /* clean up buffers */
569 mtx_zero_real64(rvec,len);
570 mtx_zero_int32(ivec,len);
571 mtx_zero_char(mark,len);
572 mtx_null_mark_release();
573 mtx_null_sum_release();
574 mtx_null_index_release();
575 /* quit */
576 return(rmax);
577 }
578
579 int32 mtx_nonzeros_in_row( mtx_matrix_t mtx, int32 row, mtx_range_t *rng)
580 {
581 struct element_t *elt;
582 int32 *tocur;
583 int32 count;
584
585 #if MTX_DEBUG
586 if(!mtx_check_matrix(mtx)) return mtx_NONE;
587 #endif
588 count = ZERO;
589 tocur = mtx->perm.col.org_to_cur;
590 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
591 if( rng == mtx_ALL_COLS )
592 for( ; NOTNULL(elt); elt = elt->next.col ) ++count;
593 else if( rng->high >= rng->low )
594 for( ; NOTNULL(elt); elt = elt->next.col )
595 if( in_range(rng,tocur[elt->col]) ) ++count;
596
597 return(count);
598 }
599
600 int32 mtx_nonzeros_in_col( mtx_matrix_t mtx, int32 col, mtx_range_t *rng)
601 {
602 struct element_t *elt;
603 int32 *tocur;
604 int32 count;
605
606 #if MTX_DEBUG
607 if(!mtx_check_matrix(mtx)) return mtx_NONE;
608 #endif
609 count = ZERO;
610 tocur = mtx->perm.row.org_to_cur;
611 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
612 if( rng == mtx_ALL_ROWS )
613 for( ; NOTNULL(elt); elt = elt->next.row ) ++count;
614 else if( rng->high >= rng->low )
615 for( ; NOTNULL(elt); elt = elt->next.row )
616 if( in_range(rng,tocur[elt->row]) ) ++count;
617
618 return(count);
619 }
620
621 int32 mtx_nonzeros_in_region(mtx_matrix_t mtx,mtx_region_t *reg)
622 {
623 struct element_t *elt;
624 int32 count;
625
626 #if MTX_DEBUG
627 if(!mtx_check_matrix(mtx)) return mtx_NONE;
628 #endif
629 count = ZERO;
630 if( reg == mtx_ENTIRE_MATRIX ) {
631 int32 org_row;
632 for( org_row = ZERO ; org_row < mtx->order ; org_row++)
633 for (elt = mtx->hdr.row[org_row] ; NOTNULL(elt) ; elt = elt->next.col)
634 count++;
635 } else {
636 if( reg->row.high >= reg->row.low && reg->col.high >= reg->col.low) {
637 int32 cur_row,rlim;
638 int32 *toorg, *tocur;
639 mtx_range_t *rng;
640 tocur = mtx->perm.col.org_to_cur;
641 toorg = mtx->perm.row.cur_to_org;
642 rlim=reg->row.high;
643 rng=&(reg->col);
644 for (cur_row=reg->row.low; cur_row<=rlim; cur_row++) {
645 elt = mtx->hdr.row[toorg[cur_row]];
646 for( ; NOTNULL(elt); elt = elt->next.col )
647 if( in_range(rng,tocur[elt->col]) ) ++count;
648 }
649 }
650 }
651 return(count);
652 }
653
654
655 int32 mtx_numbers_in_row( mtx_matrix_t mtx, int32 row, mtx_range_t *rng)
656 {
657 struct element_t *elt;
658 int32 *tocur;
659 int32 count;
660
661 #if MTX_DEBUG
662 if(!mtx_check_matrix(mtx)) return mtx_NONE;
663 #endif
664 count = ZERO;
665 tocur = mtx->perm.col.org_to_cur;
666 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
667 if( rng == mtx_ALL_COLS ) {
668 for( ; NOTNULL(elt); elt = elt->next.col ) {
669 if (elt->value != D_ZERO) ++count;
670 }
671 } else {
672 if( rng->high >= rng->low ) {
673 for( ; NOTNULL(elt); elt = elt->next.col ) {
674 if( in_range(rng,tocur[elt->col]) && elt->value!= D_ZERO ) ++count;
675 }
676 }
677 }
678 return(count);
679 }
680
681 int32 mtx_numbers_in_col( mtx_matrix_t mtx, int32 col, mtx_range_t *rng)
682 {
683 struct element_t *elt;
684 int32 *tocur;
685 int32 count;
686
687 #if MTX_DEBUG
688 if(!mtx_check_matrix(mtx)) return mtx_NONE;
689 #endif
690 count = ZERO;
691 tocur = mtx->perm.row.org_to_cur;
692 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
693 if( rng == mtx_ALL_ROWS ) {
694 for( ; NOTNULL(elt); elt = elt->next.row ) {
695 if (elt->value != D_ZERO)
696 ++count;
697 }
698 } else {
699 if( rng->high >= rng->low ) {
700 for( ; NOTNULL(elt); elt = elt->next.row ) {
701 if( in_range(rng,tocur[elt->row]) && elt->value!= D_ZERO )
702 ++count;
703 }
704 }
705 }
706 return(count);
707 }
708
709 int32 mtx_numbers_in_region(mtx_matrix_t mtx,mtx_region_t *reg)
710 {
711 struct element_t *elt;
712 int32 count;
713
714 #if MTX_DEBUG
715 if(!mtx_check_matrix(mtx)) return mtx_NONE;
716 #endif
717 count = ZERO;
718 if( reg == mtx_ENTIRE_MATRIX ) {
719 int32 org_row;
720 for( org_row = ZERO ; org_row < mtx->order ; org_row++)
721 for (elt = mtx->hdr.row[org_row] ; NOTNULL(elt) ; elt = elt->next.col)
722 if (elt->value!=D_ZERO) count++;
723 } else {
724 if( reg->row.high >= reg->row.low && reg->col.high >= reg->col.low) {
725 int32 cur_row,rlim;
726 int32 *toorg, *tocur;
727 mtx_range_t *rng;
728 tocur = mtx->perm.col.org_to_cur;
729 toorg = mtx->perm.row.cur_to_org;
730 rlim=reg->row.high;
731 rng=&(reg->col);
732 for (cur_row=reg->row.low; cur_row<=rlim; cur_row++) {
733 elt = mtx->hdr.row[toorg[cur_row]];
734 for( ; NOTNULL(elt); elt = elt->next.col )
735 if( in_range(rng,tocur[elt->col]) && elt->value!=D_ZERO ) ++count;
736 }
737 }
738 }
739 return(count);
740 }
741
742 static real64 msqr(register real64 d) {
743 return d*d;
744 }
745
746 void mtx_org_row_vec(mtx_matrix_t mtx, int32 row,
747 real64 *vec, mtx_range_t *rng)
748 {
749 struct element_t *elt;
750
751 #if MTX_DEBUG
752 if(!mtx_check_matrix(mtx)) return;
753 #endif
754 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
755 if( rng == mtx_ALL_COLS ) {
756 for( ; NOTNULL(elt); elt = elt->next.col )
757 vec[elt->col]=elt->value;
758 return;
759 }
760 if( rng->high >= rng->low ) {
761 int32 *tocur;
762 register int32 org_col, cur_col;
763 tocur = mtx->perm.col.org_to_cur;
764 for( ; NOTNULL(elt); elt = elt->next.col ) {
765 cur_col=tocur[(org_col = elt->col)];
766 if( in_range(rng,cur_col) )
767 vec[org_col]=elt->value;
768 }
769 }
770 }
771
772 void mtx_org_col_vec(mtx_matrix_t mtx, int32 col,
773 real64 *vec,mtx_range_t *rng)
774 {
775 struct element_t *elt;
776
777 #if MTX_DEBUG
778 if(!mtx_check_matrix(mtx)) return;
779 #endif
780 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
781 if( rng == mtx_ALL_ROWS ) {
782 for( ; NOTNULL(elt); elt = elt->next.row )
783 vec[elt->row]=elt->value;
784 return;
785 }
786 if( rng->high >= rng->low ) {
787 int32 *tocur;
788 register int32 org_row, cur_row;
789 tocur = mtx->perm.row.org_to_cur;
790 for( ; NOTNULL(elt); elt = elt->next.row ) {
791 cur_row=tocur[(org_row = elt->row)];
792 if( in_range(rng,cur_row) )
793 vec[org_row]=elt->value;
794 }
795 }
796 }
797
798 void mtx_cur_row_vec(mtx_matrix_t mtx, int32 row,
799 real64 *vec, mtx_range_t *rng)
800 {
801 struct element_t *elt;
802 int32 *tocur;
803
804 #if MTX_DEBUG
805 if(!mtx_check_matrix(mtx)) return;
806 #endif
807 tocur = mtx->perm.col.org_to_cur;
808 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
809 if( rng == mtx_ALL_COLS ) {
810 for( ; NOTNULL(elt); elt = elt->next.col )
811 vec[tocur[elt->col]]=elt->value;
812 return;
813 }
814 if( rng->high >= rng->low ) {
815 register int32 cur_col;
816 for( ; NOTNULL(elt); elt = elt->next.col ) {
817 cur_col=tocur[elt->col];
818 if( in_range(rng,cur_col) )
819 vec[cur_col]=elt->value;
820 }
821 }
822 }
823
824 void mtx_cur_col_vec(mtx_matrix_t mtx, int32 col,
825 real64 *vec,mtx_range_t *rng)
826 {
827 struct element_t *elt;
828 int32 *tocur;
829
830 #if MTX_DEBUG
831 if(!mtx_check_matrix(mtx)) return;
832 #endif
833 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
834 tocur = mtx->perm.row.org_to_cur;
835 if( rng == mtx_ALL_ROWS ) {
836 for( ; NOTNULL(elt); elt = elt->next.row )
837 vec[tocur[elt->row]]=elt->value;
838 return;
839 }
840 if( rng->high >= rng->low ) {
841 register int32 cur_row;
842 for( ; NOTNULL(elt); elt = elt->next.row ) {
843 cur_row=tocur[elt->row];
844 if( in_range(rng,cur_row) )
845 vec[cur_row]=elt->value;
846 }
847 }
848 }
849
850 mtx_sparse_t *mtx_org_row_sparse(mtx_matrix_t mtx, int32 row,
851 mtx_sparse_t * const vec, mtx_range_t *rng,
852 int zeroes)
853 {
854 struct element_t *elt;
855 int32 k,cap;
856
857 #if MTX_DEBUG
858 if(!mtx_check_matrix(mtx) || !mtx_check_sparse(vec)) return NULL;
859 #endif
860 if (row < 0 || row >= mtx->order) {
861 FPRINTF(g_mtxerr,"ERROR: (mtx_org_row_sparse) Row out of range.\n");
862 return NULL;
863 }
864 cap = vec->cap;
865 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
866 if( rng == mtx_ALL_COLS ) {
867 if (zeroes == mtx_IGNORE_ZEROES) {
868 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.col ) {
869 if (elt->value != D_ZERO) {
870 vec->data[k] = elt->value;
871 vec->idata[k] = elt->col;
872 k++;
873 }
874 }
875 } else {
876 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.col ) {
877 vec->data[k] = elt->value;
878 vec->idata[k] = elt->col;
879 k++;
880 }
881 }
882 vec->len = k;
883 if (ISNULL(elt)) {
884 return vec;
885 }
886 return NULL; /* ran out of vector */
887 }
888 if( rng->high >= rng->low ) {
889 int32 *tocur;
890 register int32 cur_col;
891
892 tocur = mtx->perm.col.org_to_cur;
893 if (zeroes == mtx_IGNORE_ZEROES) {
894 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.col ) {
895 if (elt->value != D_ZERO) {
896 cur_col = tocur[elt->col];
897 if( in_range(rng,cur_col) ) {
898 vec->idata[k] = elt->col;
899 vec->data[k] = elt->value;
900 k++;
901 }
902 }
903 }
904 } else {
905 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.col ) {
906 cur_col = tocur[elt->col];
907 if( in_range(rng,cur_col) ) {
908 vec->idata[k] = elt->col;
909 vec->data[k] = elt->value;
910 k++;
911 }
912 }
913 }
914 vec->len = k;
915 if (ISNULL(elt)) {
916 return vec;
917 }
918 return NULL; /* ran out of vector */
919 } else { /* nothing can exist in backward range. */
920 vec->len = 0;
921 return vec;
922 }
923 }
924
925 mtx_sparse_t *mtx_org_col_sparse(mtx_matrix_t mtx, int32 col,
926 mtx_sparse_t * const vec, mtx_range_t *rng,
927 int zeroes)
928 {
929 struct element_t *elt;
930 int32 k,cap;
931
932 #if MTX_DEBUG
933 if(!mtx_check_matrix(mtx) || !mtx_check_sparse(vec)) return NULL;
934 #endif
935 if (col < 0 || col >= mtx->order) {
936 FPRINTF(g_mtxerr,"ERROR: (mtx_org_col_sparse) Col out of range.\n");
937 return NULL;
938 }
939 cap = vec->cap;
940 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
941 if( rng == mtx_ALL_ROWS ) {
942 if (zeroes == mtx_IGNORE_ZEROES) {
943 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.row ) {
944 if (elt->value != D_ZERO) {
945 vec->data[k] = elt->value;
946 vec->idata[k] = elt->row;
947 k++;
948 }
949 }
950 } else {
951 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.row ) {
952 vec->data[k] = elt->value;
953 vec->idata[k] = elt->row;
954 k++;
955 }
956 }
957 vec->len = k;
958 if (ISNULL(elt)) {
959 return vec;
960 }
961 return NULL; /* ran out of vector */
962 }
963 if( rng->high >= rng->low ) {
964 int32 *tocur;
965 register int32 cur_row;
966
967 tocur = mtx->perm.row.org_to_cur;
968 if (zeroes == mtx_IGNORE_ZEROES) {
969 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.row ) {
970 if (elt->value != D_ZERO) {
971 cur_row = tocur[elt->row];
972 if( in_range(rng,cur_row) ) {
973 vec->idata[k] = elt->row;
974 vec->data[k] = elt->value;
975 k++;
976 }
977 }
978 }
979 } else {
980 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.row ) {
981 cur_row = tocur[elt->row];
982 if( in_range(rng,cur_row) ) {
983 vec->idata[k] = elt->row;
984 vec->data[k] = elt->value;
985 k++;
986 }
987 }
988 }
989 vec->len = k;
990 if (ISNULL(elt)) {
991 return vec;
992 }
993 return NULL; /* ran out of vector */
994 } else { /* nothing can exist in backward range. */
995 vec->len = 0;
996 return vec;
997 }
998 }
999
1000 mtx_sparse_t *mtx_cur_row_sparse(mtx_matrix_t mtx, int32 row,
1001 mtx_sparse_t * const vec, mtx_range_t *rng,
1002 int zeroes)
1003 {
1004 int32 *tocur;
1005 const struct element_t *elt;
1006 int32 cap,k;
1007
1008 #if MTX_DEBUG
1009 if(!mtx_check_matrix(mtx) || !mtx_check_sparse(vec)) return NULL;
1010 #endif
1011 if (row < 0 || row >= mtx->order) {
1012 FPRINTF(g_mtxerr,"ERROR: (mtx_cur_row_sparse) Row out of range.\n");
1013 return NULL;
1014 }
1015 tocur = mtx->perm.col.org_to_cur;
1016 cap = vec->cap;
1017 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1018 if( rng == mtx_ALL_COLS ) {
1019 if (zeroes == mtx_IGNORE_ZEROES) {
1020 for(k=0; NOTNULL(elt) && k < cap ; elt = elt->next.col ) {
1021 if (elt->value != D_ZERO) {
1022 vec->data[k] = elt->value;
1023 vec->idata[k] = tocur[elt->col];
1024 k++;
1025 }
1026 }
1027 } else {
1028 for(k=0; NOTNULL(elt) && k < cap ; elt = elt->next.col ) {
1029 vec->data[k] = elt->value;
1030 vec->idata[k] = tocur[elt->col];
1031 k++;
1032 }
1033 }
1034 vec->len = k;
1035 if (ISNULL(elt)) {
1036 return vec;
1037 }
1038 return NULL; /* ran out of vector */
1039 }
1040 /* range is not ALL */
1041 if( rng->high >= rng->low ) {
1042 register int32 cur_col;
1043
1044 if (zeroes == mtx_IGNORE_ZEROES) {
1045 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.col ) {
1046 if (elt->value != D_ZERO) {
1047 cur_col = tocur[elt->col];
1048 if( in_range(rng,cur_col) ) {
1049 vec->idata[k] = cur_col;
1050 vec->data[k] = elt->value;
1051 k++;
1052 }
1053 }
1054 }
1055 } else {
1056 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.col ) {
1057 cur_col = tocur[elt->col];
1058 if( in_range(rng,cur_col) ) {
1059 vec->idata[k] = cur_col;
1060 vec->data[k] = elt->value;
1061 k++;
1062 }
1063 }
1064 }
1065 vec->len = k;
1066 if (ISNULL(elt)) {
1067 return vec;
1068 }
1069 return NULL; /* ran out of vector */
1070 } else { /* nothing can exist in backward range. set len if needed. */
1071 vec->len = 0;
1072 return vec;
1073 }
1074 }
1075
1076 mtx_sparse_t *mtx_cur_col_sparse(mtx_matrix_t mtx, int32 col,
1077 mtx_sparse_t * const vec,mtx_range_t *rng,
1078 int zeroes)
1079 {
1080 int32 *tocur;
1081 const struct element_t *elt;
1082 int32 cap,k;
1083
1084 #if MTX_DEBUG
1085 if(!mtx_check_matrix(mtx) || !mtx_check_sparse(vec)) return NULL;
1086 #endif
1087 if (col < 0 || col >= mtx->order) {
1088 FPRINTF(g_mtxerr,"ERROR: (mtx_cur_col_sparse) Col out of range.\n");
1089 return NULL;
1090 }
1091 tocur = mtx->perm.row.org_to_cur;
1092 cap = vec->cap;
1093 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1094 if( rng == mtx_ALL_ROWS ) {
1095 if (zeroes == mtx_IGNORE_ZEROES) {
1096 for(k=0; NOTNULL(elt); elt = elt->next.row ) {
1097 if (elt->value != D_ZERO) {
1098 vec->data[k] = elt->value;
1099 vec->idata[k] = tocur[elt->row];
1100 k++;
1101 }
1102 }
1103 } else {
1104 for(k=0; NOTNULL(elt); elt = elt->next.row ) {
1105 vec->data[k] = elt->value;
1106 vec->idata[k] = tocur[elt->row];
1107 k++;
1108 }
1109 }
1110 vec->len = k;
1111 if (ISNULL(elt)) {
1112 return vec;
1113 }
1114 return NULL; /* ran out of vector */
1115 }
1116 if( rng->high >= rng->low ) {
1117 register int32 cur_row;
1118
1119 if (zeroes == mtx_IGNORE_ZEROES) {
1120 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.row ) {
1121 if (elt->value != D_ZERO) {
1122 cur_row = tocur[elt->row];
1123 if( in_range(rng,cur_row) ) {
1124 vec->idata[k] = cur_row;
1125 vec->data[k] = elt->value;
1126 k++;
1127 }
1128 }
1129 }
1130 } else {
1131 for(k=0; NOTNULL(elt) && k < cap; elt = elt->next.row ) {
1132 cur_row = tocur[elt->row];
1133 if( in_range(rng,cur_row) ) {
1134 vec->idata[k] = cur_row;
1135 vec->data[k] = elt->value;
1136 k++;
1137 }
1138 }
1139 }
1140 vec->len = k;
1141 if (ISNULL(elt)) {
1142 return vec;
1143 }
1144 return NULL; /* ran out of vector */
1145 } else { /* nothing can exist in backward range. set len if needed. */
1146 vec->len = 0;
1147 return vec;
1148 }
1149 }
1150
1151 void mtx_zr_org_vec_using_row(mtx_matrix_t mtx,int32 row,
1152 real64 *vec,mtx_range_t *rng)
1153 {
1154 struct element_t *elt;
1155
1156 #if MTX_DEBUG
1157 if(!mtx_check_matrix(mtx)) return;
1158 #endif
1159 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1160 if( rng == mtx_ALL_COLS ) {
1161 for( ; NOTNULL(elt); elt = elt->next.col )
1162 vec[elt->col]=D_ZERO;
1163 return;
1164 }
1165 if( rng->high >= rng->low ) {
1166 register int32 org_col, cur_col;
1167 int32 *tocur;
1168 tocur = mtx->perm.col.org_to_cur;
1169 for( ; NOTNULL(elt); elt = elt->next.col ) {
1170 cur_col=tocur[(org_col = elt->col)];
1171 if( in_range(rng,cur_col) )
1172 vec[org_col]=D_ZERO;
1173 }
1174 }
1175 }
1176
1177 void mtx_zr_org_vec_using_col(mtx_matrix_t mtx, int32 col,
1178 real64 *vec,mtx_range_t *rng)
1179 {
1180 struct element_t *elt;
1181
1182 #if MTX_DEBUG
1183 if(!mtx_check_matrix(mtx)) return;
1184 #endif
1185 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1186 if( rng == mtx_ALL_ROWS ) {
1187 for( ; NOTNULL(elt); elt = elt->next.row )
1188 vec[elt->row]=D_ZERO;
1189 return;
1190 }
1191 if( rng->high >= rng->low ) {
1192 register int32 org_row, cur_row;
1193 int32 *tocur;
1194 tocur = mtx->perm.row.org_to_cur;
1195 for( ; NOTNULL(elt); elt = elt->next.row ) {
1196 cur_row=tocur[(org_row = elt->row)];
1197 if( in_range(rng,cur_row) )
1198 vec[org_row]=D_ZERO;
1199 }
1200 }
1201 }
1202
1203 void mtx_zr_cur_vec_using_row(mtx_matrix_t mtx,int32 row,
1204 real64 *vec,mtx_range_t *rng)
1205 {
1206 struct element_t *elt;
1207 int32 *tocur;
1208
1209 #if MTX_DEBUG
1210 if(!mtx_check_matrix(mtx)) return;
1211 #endif
1212 tocur = mtx->perm.col.org_to_cur;
1213 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1214 if( rng == mtx_ALL_COLS ) {
1215 for( ; NOTNULL(elt); elt = elt->next.col )
1216 vec[tocur[elt->col]]=D_ZERO;
1217 return;
1218 }
1219 if( rng->high >= rng->low ) {
1220 register int32 cur_col;
1221 for( ; NOTNULL(elt); elt = elt->next.col ) {
1222 cur_col=tocur[elt->col];
1223 if( in_range(rng,cur_col) )
1224 vec[cur_col]=D_ZERO;
1225 }
1226 }
1227 }
1228
1229 void mtx_zr_cur_vec_using_col(mtx_matrix_t mtx, int32 col,
1230 real64 *vec,mtx_range_t *rng)
1231 {
1232 struct element_t *elt;
1233 int32 *tocur;
1234
1235 #if MTX_DEBUG
1236 if(!mtx_check_matrix(mtx)) return;
1237 #endif
1238 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1239 tocur = mtx->perm.row.org_to_cur;
1240 if( rng == mtx_ALL_ROWS ) {
1241 for( ; NOTNULL(elt); elt = elt->next.row )
1242 vec[tocur[elt->row]]=D_ZERO;
1243 return;
1244 }
1245 if( rng->high >= rng->low ) {
1246 register int32 cur_row;
1247 for( ; NOTNULL(elt); elt = elt->next.row ) {
1248 cur_row=tocur[elt->row];
1249 if( in_range(rng,cur_row) )
1250 vec[cur_row]=D_ZERO;
1251 }
1252 }
1253 }
1254
1255 real64 mtx_sum_sqrs_in_row( mtx_matrix_t mtx, int32 row, const mtx_range_t *rng)
1256 {
1257 struct element_t *elt;
1258 int32 *tocur;
1259 register real64 sum;
1260
1261 #if MTX_DEBUG
1262 if(!mtx_check_matrix(mtx)) return D_ZERO;
1263 #endif
1264 sum = D_ZERO;
1265 tocur = mtx->perm.col.org_to_cur;
1266 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1267 if( rng == mtx_ALL_COLS )
1268 for( ; NOTNULL(elt); elt = elt->next.col ) sum += msqr(elt->value);
1269 else if( rng->high >= rng->low )
1270 for( ; NOTNULL(elt); elt = elt->next.col )
1271 if( in_range(rng,tocur[elt->col]) ) sum += msqr(elt->value);
1272
1273 return(sum);
1274 }
1275
1276 real64 mtx_sum_sqrs_in_col( mtx_matrix_t mtx, int32 col, const mtx_range_t *rng)
1277 {
1278 struct element_t *elt;
1279 int32 *tocur;
1280 register real64 sum;
1281
1282 #if MTX_DEBUG
1283 if(!mtx_check_matrix(mtx)) return D_ZERO;
1284 #endif
1285 sum = D_ZERO;
1286 tocur = mtx->perm.row.org_to_cur;
1287 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1288 if( rng == mtx_ALL_ROWS )
1289 for( ; NOTNULL(elt); elt = elt->next.row ) sum += msqr(elt->value);
1290 else if( rng->high >= rng->low )
1291 for( ; NOTNULL(elt); elt = elt->next.row )
1292 if( in_range(rng,tocur[elt->row]) ) sum += msqr(elt->value);
1293
1294 return(sum);
1295 }
1296
1297 real64 mtx_sum_abs_in_row( mtx_matrix_t mtx, int32 row, const mtx_range_t *rng)
1298 {
1299 struct element_t *elt;
1300 int32 *tocur;
1301 register real64 sum;
1302
1303 #if MTX_DEBUG
1304 if(!mtx_check_matrix(mtx)) return D_ZERO;
1305 #endif
1306 sum = D_ZERO;
1307 tocur = mtx->perm.col.org_to_cur;
1308 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1309 if( rng == mtx_ALL_COLS )
1310 for( ; NOTNULL(elt); elt = elt->next.col ) sum+=fabs(elt->value);
1311 else if( rng->high >= rng->low )
1312 for( ; NOTNULL(elt); elt = elt->next.col )
1313 if( in_range(rng,tocur[elt->col]) ) sum+=fabs(elt->value);
1314
1315 return(sum);
1316 }
1317
1318 real64 mtx_sum_abs_in_col( mtx_matrix_t mtx, int32 col, const mtx_range_t *rng)
1319 {
1320 struct element_t *elt;
1321 int32 *tocur;
1322 register real64 sum;
1323
1324 #if MTX_DEBUG
1325 if(!mtx_check_matrix(mtx)) return D_ZERO;
1326 #endif
1327 sum = D_ZERO;
1328 tocur = mtx->perm.row.org_to_cur;
1329 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1330 if( rng == mtx_ALL_ROWS )
1331 for( ; NOTNULL(elt); elt = elt->next.row ) sum+=fabs(elt->value);
1332 else if( rng->high >= rng->low )
1333 for( ; NOTNULL(elt); elt = elt->next.row )
1334 if( in_range(rng,tocur[elt->row]) ) sum+=fabs(elt->value);
1335
1336 return(sum);
1337 }
1338
1339 real64 mtx_row_dot_full_org_vec(mtx_matrix_t mtx,
1340 int32 row,
1341 real64 *orgvec,
1342 mtx_range_t *rng,
1343 boolean transpose)
1344 {
1345 struct element_t *elt;
1346 int32 *tocur;
1347 register real64 sum;
1348
1349 #if MTX_DEBUG
1350 if(!mtx_check_matrix(mtx)) return D_ZERO;
1351 #endif
1352 sum = D_ZERO;
1353 tocur = mtx->perm.col.org_to_cur;
1354 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1355
1356 if (transpose) {
1357 int32 *toorg;
1358 toorg = mtx->perm.row.cur_to_org;
1359 if( rng == mtx_ALL_COLS ) {
1360 for( ; NOTNULL(elt); elt = elt->next.col ) {
1361 sum+=elt->value * orgvec[toorg[tocur[elt->col]]];
1362 }
1363 } else {
1364 if( rng->high >= rng->low ) {
1365 register int32 cur;
1366 for( ; NOTNULL(elt); elt = elt->next.col ) {
1367 cur=tocur[elt->col];
1368 if( in_range(rng,cur) ) {
1369 sum+=elt->value * orgvec[toorg[cur]];
1370 }
1371 }
1372 }
1373 }
1374 } else {
1375 if( rng == mtx_ALL_COLS ) {
1376 for( ; NOTNULL(elt); elt = elt->next.col ) {
1377 sum+=elt->value * orgvec[elt->col];
1378 }
1379 } else {
1380 if( rng->high >= rng->low ) {
1381 register int32 cur,org;
1382 for( ; NOTNULL(elt); elt = elt->next.col ) {
1383 cur=tocur[org=elt->col];
1384 if( in_range(rng,cur) ) {
1385 sum+=elt->value * orgvec[org];
1386 }
1387 }
1388 }
1389 }
1390 }
1391 return(sum);
1392 }
1393
1394 real64 mtx_col_dot_full_org_vec(mtx_matrix_t mtx,
1395 int32 col,
1396 real64 *orgvec,
1397 mtx_range_t *rng,
1398 boolean transpose)
1399 {
1400 struct element_t *elt;
1401 int32 *tocur;
1402 register real64 sum;
1403
1404 #if MTX_DEBUG
1405 if(!mtx_check_matrix(mtx)) return D_ZERO;
1406 #endif
1407 sum = D_ZERO;
1408 tocur = mtx->perm.row.org_to_cur;
1409 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1410
1411 if (transpose) {
1412 int32 *toorg;
1413 toorg = mtx->perm.col.cur_to_org;
1414 if( rng == mtx_ALL_ROWS ) {
1415 for( ; NOTNULL(elt); elt = elt->next.row ) {
1416 sum+=elt->value * orgvec[toorg[tocur[elt->row]]];
1417 }
1418 } else {
1419 if( rng->high >= rng->low ) {
1420 register int32 cur;
1421 for( ; NOTNULL(elt); elt = elt->next.row ) {
1422 cur=tocur[elt->row];
1423 if( in_range(rng,cur) ) {
1424 sum+=elt->value * orgvec[toorg[cur]];
1425 }
1426 }
1427 }
1428 }
1429 } else {
1430 if( rng == mtx_ALL_ROWS ) {
1431 for( ; NOTNULL(elt); elt = elt->next.row ) {
1432 sum+=elt->value * orgvec[elt->row];
1433 }
1434 } else {
1435 if( rng->high >= rng->low ) {
1436 register int32 org,cur;
1437 for( ; NOTNULL(elt); elt = elt->next.row ) {
1438 cur=tocur[org=elt->row];
1439 if( in_range(rng,cur) ) {
1440 sum+=elt->value * orgvec[org];
1441 }
1442 }
1443 }
1444 }
1445 }
1446 return(sum);
1447 }
1448
1449 real64 mtx_row_dot_full_cur_vec(mtx_matrix_t mtx,
1450 int32 row,
1451 real64 *curvec,
1452 mtx_range_t *rng,
1453 boolean transpose)
1454 {
1455 struct element_t *elt;
1456 int32 *tocur;
1457 register real64 sum;
1458
1459 #if MTX_DEBUG
1460 if(!mtx_check_matrix(mtx)) return D_ZERO;
1461 #endif
1462 sum = D_ZERO;
1463 if (!transpose) {
1464 tocur = mtx->perm.col.org_to_cur;
1465 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1466 if( rng == mtx_ALL_COLS ) {
1467 for( ; NOTNULL(elt); elt = elt->next.col ) {
1468 sum+=elt->value * curvec[tocur[elt->col]];
1469 }
1470 } else {
1471 if( rng->high >= rng->low ) {
1472 register int32 cur;
1473 for( ; NOTNULL(elt); elt = elt->next.col ) {
1474 cur=tocur[elt->col];
1475 if( in_range(rng,cur) ) {
1476 sum+=elt->value * curvec[cur];
1477 }
1478 }
1479 }
1480 }
1481 } else {
1482 FPRINTF(g_mtxerr,"(mtx) mtx_row_dot_full_cur_vec: Transpose version not\n");
1483 FPRINTF(g_mtxerr,"(mtx) done. Returning 0.0.\n");
1484 }
1485 return(sum);
1486 }
1487
1488
1489 real64 mtx_col_dot_full_cur_vec(mtx_matrix_t mtx,
1490 int32 col,
1491 real64 *curvec,
1492 mtx_range_t *rng,
1493 boolean transpose)
1494 {
1495 struct element_t *elt;
1496 int32 *tocur;
1497 register real64 sum;
1498
1499 #if MTX_DEBUG
1500 if(!mtx_check_matrix(mtx)) return D_ZERO;
1501 #endif
1502 sum = D_ZERO;
1503 if (!transpose) {
1504 tocur = mtx->perm.row.org_to_cur;
1505 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1506 if( rng == mtx_ALL_ROWS ) {
1507 for( ; NOTNULL(elt); elt = elt->next.row ) {
1508 sum+=elt->value * curvec[tocur[elt->row]];
1509 }
1510 } else {
1511 if( rng->high >= rng->low ) {
1512 register int32 cur;
1513 for( ; NOTNULL(elt); elt = elt->next.row ) {
1514 cur=tocur[elt->row];
1515 if( in_range(rng,cur) ) {
1516 sum+=elt->value * curvec[cur];
1517 }
1518 }
1519 }
1520 }
1521 } else {
1522 FPRINTF(g_mtxerr,"(mtx) mtx_col_dot_full_cur_vec: Transpose version not\n");
1523 FPRINTF(g_mtxerr,"(mtx) done. Returning 0.0.\n");
1524 }
1525 return sum;
1526 }
1527
1528 real64 mtx_row_dot_full_org_custom_vec(mtx_matrix_t mtx,
1529 mtx_matrix_t mtx2,
1530 int32 row,
1531 real64 *orgvec,
1532 mtx_range_t *rng,
1533 boolean transpose)
1534 {
1535 struct element_t *elt;
1536 int32 *tocur;
1537 register real64 sum;
1538
1539 #if MTX_DEBUG
1540 if(!mtx_check_matrix(mtx)) return D_ZERO;
1541 #endif
1542 sum = D_ZERO;
1543 tocur = mtx2->perm.col.org_to_cur;
1544 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1545
1546 if (transpose) {
1547 int32 *toorg;
1548 toorg = mtx2->perm.row.cur_to_org;
1549 if( rng == mtx_ALL_COLS ) {
1550 for( ; NOTNULL(elt); elt = elt->next.col ) {
1551 sum+=elt->value * orgvec[toorg[tocur[elt->col]]];
1552 }
1553 } else {
1554 if( rng->high >= rng->low ) {
1555 register int32 cur;
1556 for( ; NOTNULL(elt); elt = elt->next.col ) {
1557 cur=tocur[elt->col];
1558 if( in_range(rng,cur) ) {
1559 sum+=elt->value * orgvec[toorg[cur]];
1560 }
1561 }
1562 }
1563 }
1564 } else {
1565 int32 *toorg;
1566 toorg = mtx2->perm.col.cur_to_org;
1567 if( rng == mtx_ALL_COLS ) {
1568 for( ; NOTNULL(elt); elt = elt->next.col ) {
1569 sum+=elt->value * orgvec[toorg[tocur[elt->col]]];
1570 }
1571 } else {
1572 if( rng->high >= rng->low ) {
1573 register int32 cur;
1574 for( ; NOTNULL(elt); elt = elt->next.col ) {
1575 cur=tocur[elt->col];
1576 if( in_range(rng,cur) ) {
1577 sum+=elt->value * orgvec[toorg[cur]];
1578 }
1579 }
1580 }
1581 }
1582 }
1583 return(sum);
1584 }
1585
1586 real64 mtx_col_dot_full_org_custom_vec(mtx_matrix_t mtx,
1587 mtx_matrix_t mtx2,
1588 int32 col,
1589 real64 *orgvec,
1590 mtx_range_t *rng,
1591 boolean transpose)
1592 {
1593 struct element_t *elt;
1594 int32 *tocur;
1595 register real64 sum;
1596
1597 #if MTX_DEBUG
1598 if(!mtx_check_matrix(mtx)) return D_ZERO;
1599 #endif
1600 sum = D_ZERO;
1601 tocur = mtx2->perm.row.org_to_cur;
1602 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1603
1604 if (transpose) {
1605 int32 *toorg;
1606 toorg = mtx2->perm.col.cur_to_org;
1607 if( rng == mtx_ALL_ROWS ) {
1608 for( ; NOTNULL(elt); elt = elt->next.row ) {
1609 sum+=elt->value * orgvec[toorg[tocur[elt->row]]];
1610 }
1611 } else {
1612 if( rng->high >= rng->low ) {
1613 register int32 cur;
1614 for( ; NOTNULL(elt); elt = elt->next.row ) {
1615 cur=tocur[elt->row];
1616 if( in_range(rng,cur) ) {
1617 sum+=elt->value * orgvec[toorg[cur]];
1618 }
1619 }
1620 }
1621 }
1622 } else {
1623 int32 *toorg;
1624 toorg = mtx2->perm.row.cur_to_org;
1625 if( rng == mtx_ALL_ROWS ) {
1626 for( ; NOTNULL(elt); elt = elt->next.row ) {
1627 sum+=elt->value * orgvec[toorg[tocur[elt->row]]];
1628 }
1629 } else {
1630 if( rng->high >= rng->low ) {
1631 register int32 cur;
1632 for( ; NOTNULL(elt); elt = elt->next.row ) {
1633 cur=tocur[elt->row];
1634 if( in_range(rng,cur) ) {
1635 sum+=elt->value * orgvec[toorg[cur]];
1636 }
1637 }
1638 }
1639 }
1640 }
1641 return(sum);
1642 }
1643
1644
1645 void mtx_org_vec_add_row(mtx_matrix_t mtx, real64 *vec, int32 row,
1646 real64 factor, mtx_range_t *rng, boolean tr)
1647 {
1648 struct element_t *elt;
1649
1650 #if MTX_DEBUG
1651 if(!mtx_check_matrix(mtx)) return;
1652 #endif
1653
1654 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1655 if (tr) {
1656 int32 *row2org,*org2col;
1657
1658 row2org = mtx->perm.row.cur_to_org;
1659 org2col = mtx->perm.col.org_to_cur;
1660 if( rng == mtx_ALL_COLS ) {
1661 for( ; NOTNULL(elt); elt = elt->next.col ) {
1662 vec[row2org[org2col[elt->col]]] += factor * elt->value;
1663 }
1664 return;
1665 }
1666 if( rng->high >= rng->low ) {
1667 register int32 cur_col,lo,hi;
1668
1669 lo=rng->low; hi=rng->high;
1670 for( ; NOTNULL(elt); elt = elt->next.col ) {
1671 cur_col = org2col[elt->col];
1672 if( fast_in_range(lo,hi,cur_col) ) {
1673 vec[row2org[cur_col]] += factor * elt->value;
1674 }
1675 }
1676 }
1677 } else {
1678 if( rng == mtx_ALL_COLS ) {
1679 for( ; NOTNULL(elt); elt = elt->next.col ) {
1680 vec[elt->col] += factor * elt->value;
1681 }
1682 return;
1683 }
1684 if( rng->high >= rng->low ) {
1685 int32 *tocur;
1686 register int32 cur_col, org_col,lo,hi;
1687 lo=rng->low; hi=rng->high;
1688 tocur = mtx->perm.col.org_to_cur;
1689 for( ; NOTNULL(elt); elt = elt->next.col ) {
1690 cur_col=tocur[(org_col=elt->col)];
1691 if( fast_in_range(lo,hi,cur_col) ) {
1692 vec[org_col] += factor * elt->value;
1693 }
1694 }
1695 }
1696 }
1697 }
1698
1699 void mtx_org_vec_add_col(mtx_matrix_t mtx, real64 *vec, int32 col,
1700 real64 factor, mtx_range_t *rng, boolean tr)
1701 {
1702 struct element_t *elt;
1703
1704 #if MTX_DEBUG
1705 if(!mtx_check_matrix(mtx)) return;
1706 #endif
1707
1708 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1709 if (tr) {
1710 int32 *col2org,*org2row;
1711
1712 org2row = mtx->perm.row.org_to_cur;
1713 col2org = mtx->perm.col.cur_to_org;
1714 if( rng == mtx_ALL_ROWS ) {
1715 for( ; NOTNULL(elt); elt = elt->next.row ) {
1716 vec[col2org[org2row[elt->row]]] += factor * elt->value;
1717 }
1718 return;
1719 }
1720 if( rng->high >= rng->low ) {
1721 register int32 cur_row,lo,hi;
1722 lo=rng->low; hi=rng->high;
1723 for( ; NOTNULL(elt); elt = elt->next.row ) {
1724 cur_row = org2row[elt->row];
1725 if( fast_in_range(lo,hi,cur_row) ) {
1726 vec[col2org[cur_row]] += factor * elt->value;
1727 }
1728 }
1729 }
1730 } else {
1731 if( rng == mtx_ALL_ROWS ) {
1732 for( ; NOTNULL(elt); elt = elt->next.row ) {
1733 vec[elt->row] += factor * elt->value;
1734 }
1735 return;
1736 }
1737 if( rng->high >= rng->low ) {
1738 register int32 cur_row, org_row,lo,hi;
1739 int32 *tocur;
1740
1741 lo=rng->low; hi=rng->high;
1742 tocur = mtx->perm.row.org_to_cur;
1743 for( ; NOTNULL(elt); elt = elt->next.row ) {
1744 cur_row=tocur[(org_row=elt->row)];
1745 if( fast_in_range(lo,hi,cur_row) ) {
1746 vec[org_row] += factor * elt->value;
1747 }
1748 }
1749 }
1750 }
1751 }
1752
1753 void mtx_cur_vec_add_row(mtx_matrix_t mtx, real64 *vec, int32 row,
1754 real64 factor, mtx_range_t *rng, boolean tr)
1755 {
1756 struct element_t *elt;
1757 int32 *tocur;
1758
1759 #if MTX_DEBUG
1760 if(!mtx_check_matrix(mtx)) return;
1761 #endif
1762 tocur = mtx->perm.col.org_to_cur;
1763 elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
1764 if (!tr) {
1765 if( rng == mtx_ALL_COLS ) {
1766 for( ; NOTNULL(elt); elt = elt->next.col )
1767 vec[tocur[elt->col]] += factor * elt->value;
1768 return;
1769 }
1770 if( rng->high >= rng->low ) {
1771 register int32 cur_col,lo,hi;
1772 lo= rng->low; hi=rng->high;
1773 for( ; NOTNULL(elt); elt = elt->next.col ) {
1774 cur_col=tocur[elt->col];
1775 if( fast_in_range(lo,hi,cur_col) )
1776 vec[cur_col] += factor * elt->value;
1777 }
1778 }
1779 } else {
1780 FPRINTF(g_mtxerr,"(mtx) mtx_cur_vec_add_row: transpose==TRUE not done\n");
1781 /* don't know what it would mean in terms of permutations. */
1782 }
1783 }
1784
1785 void mtx_cur_vec_add_col(mtx_matrix_t mtx, real64 *vec, int32 col,
1786 real64 factor, mtx_range_t *rng, boolean tr)
1787 {
1788 struct element_t *elt;
1789 int32 *tocur;
1790
1791 #if MTX_DEBUG
1792 if(!mtx_check_matrix(mtx)) return;
1793 #endif
1794 if (!tr) {
1795 tocur = mtx->perm.row.org_to_cur;
1796 elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
1797 if( rng == mtx_ALL_ROWS ) {
1798 for( ; NOTNULL(elt); elt = elt->next.row )
1799 vec[tocur[elt->row]] += factor * elt->value;
1800 return;
1801 }
1802 if( rng->high >= rng->low ) {
1803 register int32 cur_row,lo,hi;
1804 lo= rng->low; hi=rng->high;
1805 for( ; NOTNULL(elt); elt = elt->next.row ) {
1806 cur_row=tocur[elt->row];
1807 if( fast_in_range(lo,hi,cur_row) )
1808 vec[cur_row] += factor * elt->value;
1809 }
1810 }
1811 } else {
1812 FPRINTF(g_mtxerr,"(mtx) mtx_cur_vec_add_col: transpose==TRUE not done\n");
1813 /* don't know what it would mean in terms of permutations. */
1814 }
1815 }
1816
1817 #undef __MTX_C_SEEN__

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22