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

Annotation of /trunk/ascend4/solver/mtx_linal.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, 9 months ago) by aw0a
File MIME type: text/x-csrc
File size: 15353 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.8 $
6     * Version control file: $RCSfile: mtx_linal.c,v $
7     * Date last modified: $Date: 1997/07/28 20:53:03 $
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     /*
44     * at present these assume the caller knows how long the vectors
45     * value and index must be ahead of time, or that they are of size
46     * mtx_order, which can get expensive. This may be a good functionality,
47     * but it needs some examples of why it should be supported before
48     * these are exported properly.
49     */
50     /*********************************************************************\
51     Start of New functions added by kaa. KAA_DEBUG
52     \*********************************************************************/
53     /*
54     * Export these if they work.
55     */
56     /* they don't, or more to the point, the duplicate the sparse_t functions
57     * and are a total waste of bytes. */
58     #ifdef DELETE_THIS_DEAD_FUNCTION
59     static int32 mtx_cur_col_to_csr(mtx_matrix_t mtx,
60     int32 col,
61     real64 *value,
62     int32 *index)
63     {
64     struct element_t *elt;
65     int32 *tocur;
66     int i = 0;
67    
68     elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
69     tocur = mtx->perm.row.org_to_cur;
70    
71     for( ; NOTNULL(elt); elt = elt->next.row ) {
72     index[i] = tocur[elt->row];
73     value[i] = elt->value;
74     i++;
75     }
76     return i;
77     }
78     #endif /* DELETE_THIS_DEAD_FUNCTION */
79    
80     #ifdef DELETE_THIS_DEAD_FUNCTION
81     static
82     int32 mtx_cur_col_to_csr_slow(mtx_matrix_t mtx,
83     int32 col,
84     mtx_range_t *rng,
85     real64 *value,
86     int32 *index)
87     {
88     struct element_t *elt;
89     int32 *tocur;
90     int32 cur;
91     int i = 0;
92    
93     elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
94     tocur = mtx->perm.row.org_to_cur;
95    
96     for( ; NOTNULL(elt); elt = elt->next.row ) {
97     cur = tocur[elt->row];
98     if (in_range(rng,cur)) {
99     index[i] = cur;
100     value[i] = elt->value;
101     i++;
102     }
103     }
104     return i;
105     }
106     #endif /* DELETE_THIS_DEAD_FUNCTION */
107    
108     #ifdef DELETE_THIS_DEAD_FUNCTION
109     static
110     int32 mtx_cur_row_to_csr(mtx_matrix_t mtx,
111     int32 row,
112     real64 *value,
113     int32 *index)
114     {
115     struct element_t *elt;
116     int32 *tocur;
117     int i = 0;
118    
119     tocur = mtx->perm.col.org_to_cur;
120     elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
121    
122     for( ; NOTNULL(elt); elt = elt->next.col) {
123     index[i] = tocur[elt->col];
124     value[i] = elt->value;
125     i++;
126     }
127     return i;
128     }
129     #endif /* DELETE_THIS_DEAD_FUNCTION */
130    
131     #ifdef DELETE_THIS_DEAD_FUNCTION
132     static
133     int32 mtx_cur_row_to_csr_slow(mtx_matrix_t mtx,
134     int32 row,
135     mtx_range_t *rng,
136     real64 *value,
137     int32 *index)
138     {
139     struct element_t *elt;
140     int32 *tocur;
141     int32 cur;
142     int i = 0;
143    
144     tocur = mtx->perm.col.org_to_cur;
145     elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
146    
147     for( ; NOTNULL(elt); elt = elt->next.col) {
148     cur = tocur[elt->col];
149     if (in_range(rng,cur)) {
150     index[i] = cur;
151     value[i] = elt->value;
152     i++;
153     }
154     }
155     return i;
156     }
157     #endif /* DELETE_THIS_DEAD_FUNCTION */
158    
159     /*********************************************************************\
160     End of New functions added by kaa. KAA_DEBUG
161     \*********************************************************************/
162    
163    
164     /*
165     * mtx_householder_transform_region(mtx,coef,sp,reg,droptol,transpose);
166     * add the dot product coef * u dot Transpose[u] dot A to A.
167     * real use: Anew = (I + k u.Transpose[u]) dot A
168     * aka A += u dot Transpose[u] dot A.
169     * bugs. assumes no incidence outside region in the rows of u.
170     * or its a feature: buys a LOT of speed.
171     */
172     #undef DUMELT
173     #define DUMELT ((struct element_t *)8)
174     void mtx_householder_transform_region(mtx_matrix_t mtx,
175     const real64 coef,
176     const mtx_sparse_t *sp,
177     const mtx_region_t *reg,
178     real64 droptol,
179     boolean transpose)
180     {
181     int32 *ctoorg;
182     struct element_t *elt;
183     real64 u,mult;
184     int32 ku, kr, numcols, kcol;
185     register int32 org;
186    
187     /* the following are reuseable buffers we must zero before releasing */
188     char *mv; /* mark buffer, should we need it, droptol */
189     struct element_t **ev; /* elt buffer for row addition */
190     real64 *hrow; /* u^T . A */
191     int32 *irow; /* list of nonzeros in hrow */
192    
193     #if MTX_DEBUG
194     if(!mtx_check_matrix(mtx)) return;
195     #endif
196    
197     if (coef==D_ZERO) return; /* adding 0 rather silly */
198     if (reg==mtx_ENTIRE_MATRIX ||
199     reg->col.low > reg->col.high || reg->row.low > reg->row.high) {
200     FPRINTF(g_mtxerr,"Error: (mtx.c) Bogus Householder region given.\n");
201     return;
202     }
203     if (sp==NULL) {
204     FPRINTF(g_mtxerr,"Error: (mtx.c) Bogus Householder u given.\n");
205     return;
206     }
207     if (sp->len == 0) return; /* I - k00t = I */
208     if (ISNULL(sp->idata) || ISNULL(sp->data)) {
209     FPRINTF(g_mtxerr,"Error: (mtx.c) Bogus Householder u data given.\n");
210     return;
211     }
212     if (droptol != D_ZERO || transpose) {
213     FPRINTF(g_mtxerr,
214     "Error: (mtx.c) Householder droptol and transpose not done.\n");
215     return;
216     }
217     mv = mtx_null_mark(mtx->order);
218     if (ISNULL(mv)) {
219     FPRINTF(g_mtxerr,"Error: (mtx.c) Householder. Insufficient memory.\n");
220     return;
221     }
222     ev = mtx_null_vector(mtx->order);
223     if (ISNULL(ev)) {
224     FPRINTF(g_mtxerr,"Error: (mtx.c) Householder. Insufficient memory.\n");
225     mtx_null_mark_release();
226     return;
227     }
228     hrow = mtx_null_sum(mtx->order);
229     if (ISNULL(hrow)) {
230     FPRINTF(g_mtxerr,"Error: (mtx.c) Householder. Insufficient memory.\n");
231     mtx_null_mark_release();
232     mtx_null_vector_release();
233     return;
234     }
235     irow = mtx_null_index(reg->col.high - reg->col.low +1);
236     if (ISNULL(irow)) {
237     FPRINTF(g_mtxerr,"Error: (mtx.c) Householder. Insufficient memory.\n");
238     mtx_null_mark_release();
239     mtx_null_vector_release();
240     mtx_null_sum_release();
241     return;
242     }
243     /* are we happy yet? ! */
244    
245     ctoorg = mtx->perm.col.cur_to_org;
246    
247     /* accumulate the dot products for the stuff in the region.no range check! */
248     /* in I- tau u u^T . A this is hrow=u^T.A. idata[ku] is an org row index. */
249     for (ku=0; ku < sp->len; ku++) {
250     u = sp->data[ku];
251     if (u != D_ZERO) {
252     elt = mtx->hdr.row[sp->idata[ku]];
253     while (NOTNULL(elt)) {
254     hrow[elt->col] += u*elt->value; /* 15% */
255     elt = elt->next.col; /* 4% */
256     }
257     }
258     }
259     /* accumulate the indices needed to drive zeroing */
260     kr = 0;
261     for (ku = reg->col.low; ku <= reg->col.high; ku++) {
262     org = ctoorg[ku];
263     if (hrow[org]!=D_ZERO) {
264     irow[kr] = org; /* collect index */
265     ev[org] = DUMELT; /* init elt array with bogus ptr */
266     kr++;
267     }
268     }
269     /* irow 0..kr-1 now has the org col indexes of nonzero elements in hrow,ev */
270     numcols = kr;
271     /* now add hh transform to rows of A in u, cases with and without coef = 1 */
272     if (coef == D_ONE) {
273     for (ku=0; ku < sp->len; ku++) {
274     if (sp->data[ku] != D_ZERO) {
275     mult = -sp->data[ku];
276     org = sp->idata[ku];
277     /* expand row of interest into locations of interest */
278     elt = mtx->hdr.row[org];
279     while (NOTNULL(elt)) {
280     if (NOTNULL(ev[elt->col])) ev[elt->col] = elt; /* 11% */
281     elt = elt->next.col; /* 4% */
282     }
283     /* run through irow doing the math (finally) */
284     for (kr=0; kr < numcols; kr++) {
285     kcol = irow[kr]; /* 2 */
286     if (ev[kcol] > DUMELT) { /* 12% */
287     ev[kcol]->value += mult*hrow[kcol]; /* 14% */
288     ev[kcol] = DUMELT; /* 9% */
289     /* reinit ev col */
290     } else {
291     mtx_create_element_value(mtx, org, kcol, mult*hrow[kcol]); /*6*/
292     /* let ev[irow[kr]] col stay DUMELT */
293     }
294     }
295     }
296     }
297     } else {
298     for (ku=0; ku < sp->len; ku++) {
299     mult = -sp->data[ku]*coef;
300     if (mult != D_ZERO) {
301     org = sp->idata[ku];
302     /* expand row of interest into locations of interest */
303     elt = mtx->hdr.row[org];
304     while (NOTNULL(elt)) {
305     if (NOTNULL(ev[elt->col])) ev[elt->col] = elt;
306     elt = elt->next.col;
307     }
308     /* run through irow doing the math (finally) */
309     for (kr=0; kr < numcols; kr++) {
310     kcol = irow[kr];
311     if (ev[kcol] > DUMELT) {
312     ev[kcol]->value += mult*hrow[kcol];
313     ev[kcol] = DUMELT;
314     /* reinit ev col */
315     } else {
316     mtx_create_element_value(mtx, org, kcol, mult*hrow[kcol]);
317     /* let ev[irow[kr]] col stay DUMELT */
318     }
319     }
320     }
321     }
322     }
323     /* end case coef != 1 */
324     for (kr=0; kr < numcols; kr++) {
325     ev[irow[kr]] = NULL;
326     hrow[irow[kr]] = D_ZERO;
327     irow[kr] = 0;
328     }
329     /* pack up and go home */
330     mtx_null_mark_release();
331     mtx_null_vector_release();
332     mtx_null_sum_release();
333     mtx_null_index_release();
334     return;
335     }
336     #undef DUMELT
337    
338    
339     /*
340     *********************************************************************
341     * mtx_eliminate_row.
342     *
343     * Start of some new routines, so that we can do so micro
344     * management of the elimination routines.
345     *********************************************************************
346     */
347     struct mtx_linklist {
348     int32 index;
349     struct mtx_linklist *prev;
350     };
351    
352     /*
353     * This function assumes that head is not NULL.
354     * we dont check !!
355     */
356     #ifdef THIS_IS_AN_UNUSED_FUNCTION
357     static struct mtx_linklist *insert_link(struct mtx_linklist *head,
358     int k,
359     mem_store_t eltbuffer)
360     {
361     struct mtx_linklist *ptr1, *ptr2;
362     struct mtx_linklist *target;
363    
364     target = (struct mtx_linklist *)mem_get_element(eltbuffer);
365     target->index = k;
366    
367     ptr2 = ptr1 = head;
368     ptr1 = head->prev;
369     while (ptr1) {
370     if (ptr1->index < k) {
371     target->prev = ptr1;
372     ptr2->prev = target;
373     return target;
374     }
375     ptr2 = ptr1;
376     ptr1 = ptr1->prev;
377     }
378     /*
379     * If we are here then we reached the end of
380     * the chain. Just add target and quit.
381     */
382     target->prev = ptr1;
383     ptr2->prev = target;
384     return target;
385     }
386     #endif /* THIS_IS_AN_UNUSED_FUNCTION */
387    
388     /*
389     * head may be NULL for this function. We simply add
390     * the new index *unsorted* and return the end of top
391     * of the chain.
392     */
393     #ifdef THIS_IS_AN_UNUSED_FUNCTION
394     static struct mtx_linklist *add_link(struct mtx_linklist *head,
395     int k,
396     mem_store_t eltbuffer)
397     {
398     struct mtx_linklist *target;
399    
400     target = (struct mtx_linklist *)mem_get_element(eltbuffer);
401     target->index = k;
402    
403     target->prev = head;
404     return target;
405     }
406     #endif /* THIS_IS_AN_UNUSED_FUNCTION */
407    
408     /*
409     * NOTE: This function needs a droptol as a parameter
410     */
411    
412     #ifdef THIS_IS_AN_UNUSED_FUNCTION
413     static
414     void mtx_eliminate_row2(mtx_matrix_t mtx,
415     mtx_matrix_t upper_mtx,
416     mtx_range_t *rng,
417     int32 row,
418     real64 *vec,
419     real64 *pivots,
420     int32 *inserted,
421     mem_store_t eltbuffer)
422     {
423     struct element_t *elt;
424     struct mtx_linklist *diag, *right=NULL, *ptr;
425     int32 *tocur;
426     real64 multiplier;
427     int k, j;
428     mtx_coord_t nz;
429    
430     (void)rng; /* stop gcc whine about unused parameter */
431    
432     /*
433     * Move all non-zeros from current row to full array.
434     * The full array would have been initialized before,
435     * we must put it back in the clean state when we leave.
436     * All operations are done over mtx_ALL_COLS.
437     */
438     /*
439     * we use the following code fragment instead of :
440     * mtx_cur_row_vec(mtx,row,vec,mtx_ALL_COLS);
441     * so that we can put the elements in the correct place.
442     */
443    
444     diag = (struct mtx_linklist *)mem_get_element(eltbuffer);
445     diag->index = row;
446     diag->prev = NULL;
447     inserted[row] = 1;
448    
449     /*
450     * the insertion for this phase.
451     */
452     tocur = mtx->perm.col.org_to_cur;
453     elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
454     ptr = diag;
455     for ( ;NOTNULL(elt); elt = elt->next.col) {
456     if (elt->value == 0.0) {
457     continue; /* hard zeros */
458     }
459     k = tocur[elt->col];
460     vec[k] = elt->value;
461     if (k < row) { /* the less than is critical */
462     (void)insert_link(diag,k,eltbuffer);
463     inserted[k] = 1;
464     }
465     else if (k > row) {
466     right = add_link(right,k,eltbuffer);
467     inserted[k] = 1;
468     }
469     else
470     continue;
471     }
472    
473     mtx_clear_row(mtx,row,mtx_ALL_COLS);
474    
475     /* we shuold be trapping for these 0 multipliers before. !!! */
476    
477     for (ptr = diag->prev; NOTNULL(ptr); ) {
478     k = ptr->index;
479     multiplier = vec[k]/pivots[k];
480     if (multiplier==D_ZERO) {
481     ptr = ptr->prev;
482     /* FPRINTF(stderr,"0 multiplier found at %d\n",k); */
483     continue;
484     }
485     #ifdef NOP_DEBUG
486     mtx_number_ops++;
487     #endif /* NOP_DEBUG */
488     elt = mtx->hdr.row[mtx->perm.row.cur_to_org[k]];
489     for ( ;NOTNULL(elt); elt = elt->next.col) {
490     j = tocur[elt->col];
491     if (!inserted[j]) {
492     if (j < k)
493     (void)insert_link(ptr,j,eltbuffer);
494     else
495     right = add_link(right,j,eltbuffer);
496     inserted[j] = 1;
497     }
498     vec[j] = vec[j] - multiplier * elt->value;
499     #ifdef NOP_DEBUG
500     mtx_number_ops++;
501     #endif /* NOP_DEBUG */
502     }
503     vec[k] = multiplier; /* backpatch multiplier */
504     ptr = ptr->prev;
505     }
506    
507     /*
508     * Map the data back to the appropriate matrices.
509     */
510    
511     /*
512     * Fill up the upper_matrix with the multipliers.
513     */
514     nz.row = row;
515     for (ptr = diag->prev; NOTNULL(ptr); ptr = ptr->prev) {
516     nz.col = ptr->index;
517     if (vec[nz.col] != D_ZERO) { /* dont fill hard 0's */
518     mtx_fill_value(upper_mtx,&nz,vec[nz.col]);
519     }
520     vec[nz.col] = D_ZERO;
521     inserted[nz.col] = 0;
522     }
523    
524     /*
525     * Fill the diagonal back to the regular matrix.
526     */
527     nz.col = row;
528     if (vec[nz.col] != D_ZERO) { /* dont fill hard 0's */
529     mtx_fill_value(mtx,&nz,vec[nz.col]);
530     }
531     vec[row] = D_ZERO;
532     inserted[row] = 0;
533    
534     /*
535     * Fill the lower matrix with the stuff to the right of
536     * diagonal.
537     */
538     for (ptr = right; NOTNULL(ptr); ptr = ptr->prev) {
539     nz.col = ptr->index;
540     if (fabs(vec[nz.col]) > 1.0e-16) { /* THIS NEEDS A DROP TOL DEFINE */
541     mtx_fill_value(mtx,&nz,vec[nz.col]);
542     }
543     vec[nz.col] = D_ZERO;
544     inserted[nz.col] = 0;
545     }
546     }
547     #endif /* THIS_IS_AN_UNUSED_FUNCTION */
548    
549     #undef __MTX_C_SEEN__

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