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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 8 months ago) by aw0a
File MIME type: text/x-csrc
File size: 50034 byte(s)
Setting up web subdirectory in repository
1 aw0a 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