home *** CD-ROM | disk | FTP | other *** search
Text File | 1994-06-07 | 176.0 KB | 7,044 lines | [TEXT/ttxt] |
- # to unbundle, sh this file (in an empty directory)
- echo sparse.c 1>&2
- sed >sparse.c <<'//GO.SYSIN DD sparse.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -/*
- - Sparse matrix package
- - See also: sparse.h, matrix.h
- - */
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include <stdlib.h>
- -#include "sparse.h"
- -
- -
- -static char rcsid[] = "$Id: sparse.c,v 1.10 1994/03/08 05:46:07 des Exp $";
- -
- -#define MINROWLEN 10
- -
- -
- -
- -/* sp_get_val -- returns the (i,j) entry of the sparse matrix A */
- -double sp_get_val(A,i,j)
- -SPMAT *A;
- -int i, j;
- -{
- - SPROW *r;
- - int idx;
- -
- - if ( A == SMNULL )
- - error(E_NULL,"sp_get_val");
- - if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
- - error(E_SIZES,"sp_get_val");
- -
- - r = A->row+i;
- - idx = sprow_idx(r,j);
- - if ( idx < 0 )
- - return 0.0;
- - /* else */
- - return r->elt[idx].val;
- -}
- -
- -/* sp_set_val -- sets the (i,j) entry of the sparse matrix A */
- -double sp_set_val(A,i,j,val)
- -SPMAT *A;
- -int i, j;
- -double val;
- -{
- - SPROW *r;
- - int idx, idx2, new_len;
- -
- - if ( A == SMNULL )
- - error(E_NULL,"sp_set_val");
- - if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
- - error(E_SIZES,"sp_set_val");
- -
- - r = A->row+i;
- - idx = sprow_idx(r,j);
- - /* printf("sp_set_val: idx = %d\n",idx); */
- - if ( idx >= 0 )
- - { r->elt[idx].val = val; return val; }
- - /* else */ if ( idx < -1 )
- - {
- - /* Note: this destroys the column & diag access paths */
- - A->flag_col = A->flag_diag = FALSE;
- - /* shift & insert new value */
- - idx = -(idx+2); /* this is the intended insertion index */
- - if ( r->len >= r->maxlen )
- - {
- - r->len = r->maxlen;
- - new_len = max(2*r->maxlen+1,5);
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,A->row[i].maxlen*sizeof(row_elt),
- - new_len*sizeof(row_elt));
- - }
- -
- - r->elt = RENEW(r->elt,new_len,row_elt);
- - if ( ! r->elt ) /* can't allocate */
- - error(E_MEM,"sp_set_val");
- - r->maxlen = 2*r->maxlen+1;
- - }
- - for ( idx2 = r->len-1; idx2 >= idx; idx2-- )
- - MEM_COPY((char *)(&(r->elt[idx2])),
- - (char *)(&(r->elt[idx2+1])),sizeof(row_elt));
- - /************************************************************
- - if ( idx < r->len )
- - MEM_COPY((char *)(&(r->elt[idx])),(char *)(&(r->elt[idx+1])),
- - (r->len-idx)*sizeof(row_elt));
- - ************************************************************/
- - r->len++;
- - r->elt[idx].col = j;
- - return r->elt[idx].val = val;
- - }
- - /* else -- idx == -1, error in index/matrix! */
- - return 0.0;
- -}
- -
- -/* sp_mv_mlt -- sparse matrix/dense vector multiply
- - -- result is in out, which is returned unless out==NULL on entry
- - -- if out==NULL on entry then the result vector is created */
- -VEC *sp_mv_mlt(A,x,out)
- -SPMAT *A;
- -VEC *x, *out;
- -{
- - int i, j_idx, m, n, max_idx;
- - Real sum, *x_ve;
- - SPROW *r;
- - row_elt *elts;
- -
- - if ( ! A || ! x )
- - error(E_NULL,"sp_mv_mlt");
- - if ( x->dim != A->n )
- - error(E_SIZES,"sp_mv_mlt");
- - if ( ! out || out->dim < A->m )
- - out = v_resize(out,A->m);
- - if ( out == x )
- - error(E_INSITU,"sp_mv_mlt");
- - m = A->m; n = A->n;
- - x_ve = x->ve;
- -
- - for ( i = 0; i < m; i++ )
- - {
- - sum = 0.0;
- - r = &(A->row[i]);
- - max_idx = r->len;
- - elts = r->elt;
- - for ( j_idx = 0; j_idx < max_idx; j_idx++, elts++ )
- - sum += elts->val*x_ve[elts->col];
- - out->ve[i] = sum;
- - }
- - return out;
- -}
- -
- -/* sp_vm_mlt -- sparse matrix/dense vector multiply from left
- - -- result is in out, which is returned unless out==NULL on entry
- - -- if out==NULL on entry then result vector is created & returned */
- -VEC *sp_vm_mlt(A,x,out)
- -SPMAT *A;
- -VEC *x, *out;
- -{
- - int i, j_idx, m, n, max_idx;
- - Real tmp, *x_ve, *out_ve;
- - SPROW *r;
- - row_elt *elts;
- -
- - if ( ! A || ! x )
- - error(E_NULL,"sp_vm_mlt");
- - if ( x->dim != A->m )
- - error(E_SIZES,"sp_vm_mlt");
- - if ( ! out || out->dim < A->n )
- - out = v_resize(out,A->n);
- - if ( out == x )
- - error(E_INSITU,"sp_vm_mlt");
- -
- - m = A->m; n = A->n;
- - v_zero(out);
- - x_ve = x->ve; out_ve = out->ve;
- -
- - for ( i = 0; i < m; i++ )
- - {
- - r = A->row+i;
- - max_idx = r->len;
- - elts = r->elt;
- - tmp = x_ve[i];
- - for ( j_idx = 0; j_idx < max_idx; j_idx++, elts++ )
- - out_ve[elts->col] += elts->val*tmp;
- - }
- -
- - return out;
- -}
- -
- -
- -/* sp_get -- get sparse matrix
- - -- len is number of elements available for each row without
- - allocating further memory */
- -SPMAT *sp_get(m,n,maxlen)
- -int m, n, maxlen;
- -{
- - SPMAT *A;
- - SPROW *rows;
- - int i;
- -
- - if ( m < 0 || n < 0 )
- - error(E_NEG,"sp_get");
- -
- - maxlen = max(maxlen,1);
- -
- - A = NEW(SPMAT);
- - if ( ! A ) /* can't allocate */
- - error(E_MEM,"sp_get");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,sizeof(SPMAT));
- - mem_numvar(TYPE_SPMAT,1);
- - }
- - /* fprintf(stderr,"Have SPMAT structure\n"); */
- -
- - A->row = rows = NEW_A(m,SPROW);
- - if ( ! A->row ) /* can't allocate */
- - error(E_MEM,"sp_get");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,m*sizeof(SPROW));
- - }
- - /* fprintf(stderr,"Have row structure array\n"); */
- -
- - A->start_row = NEW_A(n,int);
- - A->start_idx = NEW_A(n,int);
- - if ( ! A->start_row || ! A->start_idx ) /* can't allocate */
- - error(E_MEM,"sp_get");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,2*n*sizeof(int));
- - }
- - for ( i = 0; i < n; i++ )
- - A->start_row[i] = A->start_idx[i] = -1;
- - /* fprintf(stderr,"Have start_row array\n"); */
- -
- - A->m = A->max_m = m;
- - A->n = A->max_n = n;
- -
- - for ( i = 0; i < m; i++, rows++ )
- - {
- - rows->elt = NEW_A(maxlen,row_elt);
- - if ( ! rows->elt )
- - error(E_MEM,"sp_get");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,maxlen*sizeof(row_elt));
- - }
- - /* fprintf(stderr,"Have row %d element array\n",i); */
- - rows->len = 0;
- - rows->maxlen = maxlen;
- - rows->diag = -1;
- - }
- -
- - return A;
- -}
- -
- -
- -/* sp_free -- frees up the memory for a sparse matrix */
- -int sp_free(A)
- -SPMAT *A;
- -{
- - SPROW *r;
- - int i;
- -
- - if ( ! A )
- - return -1;
- - if ( A->start_row != (int *)NULL ) {
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,A->max_n*sizeof(int),0);
- - }
- - free((char *)(A->start_row));
- - }
- - if ( A->start_idx != (int *)NULL ) {
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,A->max_n*sizeof(int),0);
- - }
- -
- - free((char *)(A->start_idx));
- - }
- - if ( ! A->row )
- - {
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,sizeof(SPMAT),0);
- - mem_numvar(TYPE_SPMAT,-1);
- - }
- -
- - free((char *)A);
- - return 0;
- - }
- - for ( i = 0; i < A->m; i++ )
- - {
- - r = &(A->row[i]);
- - if ( r->elt != (row_elt *)NULL ) {
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,A->row[i].maxlen*sizeof(row_elt),0);
- - }
- - free((char *)(r->elt));
- - }
- - }
- -
- - if (mem_info_is_on()) {
- - if (A->row)
- - mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),0);
- - mem_bytes(TYPE_SPMAT,sizeof(SPMAT),0);
- - mem_numvar(TYPE_SPMAT,-1);
- - }
- -
- - free((char *)(A->row));
- - free((char *)A);
- -
- - return 0;
- -}
- -
- -
- -/* sp_copy -- constructs a copy of a given matrix
- - -- note that the max_len fields (etc) are no larger in the copy
- - than necessary
- - -- result is returned */
- -SPMAT *sp_copy(A)
- -SPMAT *A;
- -{
- - SPMAT *out;
- - SPROW *row1, *row2;
- - int i;
- -
- - if ( A == SMNULL )
- - error(E_NULL,"sp_copy");
- - if ( ! (out=NEW(SPMAT)) )
- - error(E_MEM,"sp_copy");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,sizeof(SPMAT));
- - mem_numvar(TYPE_SPMAT,1);
- - }
- - out->m = out->max_m = A->m; out->n = out->max_n = A->n;
- -
- - /* set up rows */
- - if ( ! (out->row=NEW_A(A->m,SPROW)) )
- - error(E_MEM,"sp_copy");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,A->m*sizeof(SPROW));
- - }
- - for ( i = 0; i < A->m; i++ )
- - {
- - row1 = &(A->row[i]);
- - row2 = &(out->row[i]);
- - if ( ! (row2->elt=NEW_A(max(row1->len,3),row_elt)) )
- - error(E_MEM,"sp_copy");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,max(row1->len,3)*sizeof(row_elt));
- - }
- - row2->len = row1->len;
- - row2->maxlen = max(row1->len,3);
- - row2->diag = row1->diag;
- - MEM_COPY((char *)(row1->elt),(char *)(row2->elt),
- - row1->len*sizeof(row_elt));
- - }
- -
- - /* set up start arrays -- for column access */
- - if ( ! (out->start_idx=NEW_A(A->n,int)) ||
- - ! (out->start_row=NEW_A(A->n,int)) )
- - error(E_MEM,"sp_copy");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,2*A->n*sizeof(int));
- - }
- - MEM_COPY((char *)(A->start_idx),(char *)(out->start_idx),
- - A->n*sizeof(int));
- - MEM_COPY((char *)(A->start_row),(char *)(out->start_row),
- - A->n*sizeof(int));
- -
- - return out;
- -}
- -
- -/* sp_col_access -- set column access path; i.e. nxt_row, nxt_idx fields
- - -- returns A */
- -SPMAT *sp_col_access(A)
- -SPMAT *A;
- -{
- - int i, j, j_idx, len, m, n;
- - SPROW *row;
- - row_elt *r_elt;
- - int *start_row, *start_idx;
- -
- - if ( A == SMNULL )
- - error(E_NULL,"sp_col_access");
- -
- - m = A->m; n = A->n;
- -
- - /* initialise start_row and start_idx */
- - start_row = A->start_row; start_idx = A->start_idx;
- - for ( j = 0; j < n; j++ )
- - { *start_row++ = -1; *start_idx++ = -1; }
- -
- - start_row = A->start_row; start_idx = A->start_idx;
- -
- - /* now work UP the rows, setting nxt_row, nxt_idx fields */
- - for ( i = m-1; i >= 0; i-- )
- - {
- - row = &(A->row[i]);
- - r_elt = row->elt;
- - len = row->len;
- - for ( j_idx = 0; j_idx < len; j_idx++, r_elt++ )
- - {
- - j = r_elt->col;
- - r_elt->nxt_row = start_row[j];
- - r_elt->nxt_idx = start_idx[j];
- - start_row[j] = i;
- - start_idx[j] = j_idx;
- - }
- - }
- -
- - A->flag_col = TRUE;
- - return A;
- -}
- -
- -/* sp_diag_access -- set diagonal access path(s) */
- -SPMAT *sp_diag_access(A)
- -SPMAT *A;
- -{
- - int i, m;
- - SPROW *row;
- -
- - if ( A == SMNULL )
- - error(E_NULL,"sp_diag_access");
- -
- - m = A->m;
- -
- - row = A->row;
- - for ( i = 0; i < m; i++, row++ )
- - row->diag = sprow_idx(row,i);
- -
- - A->flag_diag = TRUE;
- -
- - return A;
- -}
- -
- -/* sp_m2dense -- convert a sparse matrix to a dense one */
- -MAT *sp_m2dense(A,out)
- -SPMAT *A;
- -MAT *out;
- -{
- - int i, j_idx;
- - SPROW *row;
- - row_elt *elt;
- -
- - if ( ! A )
- - error(E_NULL,"sp_m2dense");
- - if ( ! out || out->m < A->m || out->n < A->n )
- - out = m_get(A->m,A->n);
- -
- - m_zero(out);
- - for ( i = 0; i < A->m; i++ )
- - {
- - row = &(A->row[i]);
- - elt = row->elt;
- - for ( j_idx = 0; j_idx < row->len; j_idx++, elt++ )
- - out->me[i][elt->col] = elt->val;
- - }
- -
- - return out;
- -}
- -
- -
- -/* C = A+B, can be in situ */
- -SPMAT *sp_add(A,B,C)
- -SPMAT *A, *B, *C;
- -{
- - int i, in_situ;
- - SPROW *rc;
- - static SPROW *tmp;
- -
- - if ( ! A || ! B )
- - error(E_NULL,"sp_add");
- - if ( A->m != B->m || A->n != B->n )
- - error(E_SIZES,"sp_add");
- - if (C == A || C == B)
- - in_situ = TRUE;
- - else in_situ = FALSE;
- -
- - if ( ! C )
- - C = sp_get(A->m,A->n,5);
- - else {
- - if ( C->m != A->m || C->n != A->n )
- - error(E_SIZES,"sp_add");
- - if (!in_situ) sp_zero(C);
- - }
- -
- - if (tmp == (SPROW *)NULL && in_situ) {
- - tmp = sprow_get(MINROWLEN);
- - MEM_STAT_REG(tmp,TYPE_SPROW);
- - }
- -
- - if (in_situ)
- - for (i=0; i < A->m; i++) {
- - rc = &(C->row[i]);
- - sprow_add(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
- - sprow_resize(rc,tmp->len,TYPE_SPMAT);
- - MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
- - rc->len = tmp->len;
- - }
- - else
- - for (i=0; i < A->m; i++) {
- - sprow_add(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
- - }
- -
- - C->flag_col = C->flag_diag = FALSE;
- -
- - return C;
- -}
- -
- -/* C = A-B, cannot be in situ */
- -SPMAT *sp_sub(A,B,C)
- -SPMAT *A, *B, *C;
- -{
- - int i, in_situ;
- - SPROW *rc;
- - static SPROW *tmp;
- -
- - if ( ! A || ! B )
- - error(E_NULL,"sp_sub");
- - if ( A->m != B->m || A->n != B->n )
- - error(E_SIZES,"sp_sub");
- - if (C == A || C == B)
- - in_situ = TRUE;
- - else in_situ = FALSE;
- -
- - if ( ! C )
- - C = sp_get(A->m,A->n,5);
- - else {
- - if ( C->m != A->m || C->n != A->n )
- - error(E_SIZES,"sp_sub");
- - if (!in_situ) sp_zero(C);
- - }
- -
- - if (tmp == (SPROW *)NULL && in_situ) {
- - tmp = sprow_get(MINROWLEN);
- - MEM_STAT_REG(tmp,TYPE_SPROW);
- - }
- -
- - if (in_situ)
- - for (i=0; i < A->m; i++) {
- - rc = &(C->row[i]);
- - sprow_sub(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
- - sprow_resize(rc,tmp->len,TYPE_SPMAT);
- - MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
- - rc->len = tmp->len;
- - }
- - else
- - for (i=0; i < A->m; i++) {
- - sprow_sub(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
- - }
- -
- - C->flag_col = C->flag_diag = FALSE;
- -
- - return C;
- -}
- -
- -/* C = A+alpha*B, cannot be in situ */
- -SPMAT *sp_mltadd(A,B,alpha,C)
- -SPMAT *A, *B, *C;
- -double alpha;
- -{
- - int i, in_situ;
- - SPROW *rc;
- - static SPROW *tmp;
- -
- - if ( ! A || ! B )
- - error(E_NULL,"sp_mltadd");
- - if ( A->m != B->m || A->n != B->n )
- - error(E_SIZES,"sp_mltadd");
- - if (C == A || C == B)
- - in_situ = TRUE;
- - else in_situ = FALSE;
- -
- - if ( ! C )
- - C = sp_get(A->m,A->n,5);
- - else {
- - if ( C->m != A->m || C->n != A->n )
- - error(E_SIZES,"sp_mltadd");
- - if (!in_situ) sp_zero(C);
- - }
- -
- - if (tmp == (SPROW *)NULL && in_situ) {
- - tmp = sprow_get(MINROWLEN);
- - MEM_STAT_REG(tmp,TYPE_SPROW);
- - }
- -
- - if (in_situ)
- - for (i=0; i < A->m; i++) {
- - rc = &(C->row[i]);
- - sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,tmp,TYPE_SPROW);
- - sprow_resize(rc,tmp->len,TYPE_SPMAT);
- - MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
- - rc->len = tmp->len;
- - }
- - else
- - for (i=0; i < A->m; i++) {
- - sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,
- - &(C->row[i]),TYPE_SPMAT);
- - }
- -
- - C->flag_col = C->flag_diag = FALSE;
- -
- - return C;
- -}
- -
- -
- -
- -/* B = alpha*A, can be in situ */
- -SPMAT *sp_smlt(A,alpha,B)
- -SPMAT *A, *B;
- -double alpha;
- -{
- - int i;
- -
- - if ( ! A )
- - error(E_NULL,"sp_smlt");
- - if ( ! B )
- - B = sp_get(A->m,A->n,5);
- - else
- - if ( A->m != B->m || A->n != B->n )
- - error(E_SIZES,"sp_smlt");
- -
- - for (i=0; i < A->m; i++) {
- - sprow_smlt(&(A->row[i]),alpha,0,&(B->row[i]),TYPE_SPMAT);
- - }
- - return B;
- -}
- -
- -
- -
- -/* sp_zero -- zero all the (represented) elements of a sparse matrix */
- -SPMAT *sp_zero(A)
- -SPMAT *A;
- -{
- - int i, idx, len;
- - row_elt *elt;
- -
- - if ( ! A )
- - error(E_NULL,"sp_zero");
- -
- - for ( i = 0; i < A->m; i++ )
- - {
- - elt = A->row[i].elt;
- - len = A->row[i].len;
- - for ( idx = 0; idx < len; idx++ )
- - (*elt++).val = 0.0;
- - }
- -
- - return A;
- -}
- -
- -/* sp_copy2 -- copy sparse matrix (type 2)
- - -- keeps structure of the OUT matrix */
- -SPMAT *sp_copy2(A,OUT)
- -SPMAT *A, *OUT;
- -{
- - int i /* , idx, len1, len2 */;
- - SPROW *r1, *r2;
- - static SPROW *scratch = (SPROW *)NULL;
- - /* row_elt *e1, *e2; */
- -
- - if ( ! A )
- - error(E_NULL,"sp_copy2");
- - if ( ! OUT )
- - OUT = sp_get(A->m,A->n,10);
- - if ( ! scratch ) {
- - scratch = sprow_xpd(scratch,MINROWLEN,TYPE_SPROW);
- - MEM_STAT_REG(scratch,TYPE_SPROW);
- - }
- -
- - if ( OUT->m < A->m )
- - {
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
- - A->m*sizeof(SPROW));
- - }
- -
- - OUT->row = RENEW(OUT->row,A->m,SPROW);
- - if ( ! OUT->row )
- - error(E_MEM,"sp_copy2");
- -
- - for ( i = OUT->m; i < A->m; i++ )
- - {
- - OUT->row[i].elt = NEW_A(MINROWLEN,row_elt);
- - if ( ! OUT->row[i].elt )
- - error(E_MEM,"sp_copy2");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
- - }
- -
- - OUT->row[i].maxlen = MINROWLEN;
- - OUT->row[i].len = 0;
- - }
- - OUT->m = A->m;
- - }
- -
- - OUT->flag_col = OUT->flag_diag = FALSE;
- - /* sp_zero(OUT); */
- -
- - for ( i = 0; i < A->m; i++ )
- - {
- - r1 = &(A->row[i]); r2 = &(OUT->row[i]);
- - sprow_copy(r1,r2,scratch,TYPE_SPROW);
- - if ( r2->maxlen < scratch->len )
- - sprow_xpd(r2,scratch->len,TYPE_SPMAT);
- - MEM_COPY((char *)(scratch->elt),(char *)(r2->elt),
- - scratch->len*sizeof(row_elt));
- - r2->len = scratch->len;
- - /*******************************************************
- - e1 = r1->elt; e2 = r2->elt;
- - len1 = r1->len; len2 = r2->len;
- - for ( idx = 0; idx < len2; idx++, e2++ )
- - e2->val = 0.0;
- - for ( idx = 0; idx < len1; idx++, e1++ )
- - sprow_set_val(r2,e1->col,e1->val);
- - *******************************************************/
- - }
- -
- - sp_col_access(OUT);
- - return OUT;
- -}
- -
- -/* sp_resize -- resize a sparse matrix
- - -- don't destroying any contents if possible
- - -- returns resized matrix */
- -SPMAT *sp_resize(A,m,n)
- -SPMAT *A;
- -int m, n;
- -{
- - int i, len;
- - SPROW *r;
- -
- - if (m < 0 || n < 0)
- - error(E_NEG,"sp_resize");
- -
- - if ( ! A )
- - return sp_get(m,n,10);
- -
- - if (m == A->m && n == A->n)
- - return A;
- -
- - if ( m <= A->max_m )
- - {
- - for ( i = A->m; i < m; i++ )
- - A->row[i].len = 0;
- - A->m = m;
- - }
- - else
- - {
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
- - m*sizeof(SPROW));
- - }
- -
- - A->row = RENEW(A->row,(unsigned)m,SPROW);
- - if ( ! A->row )
- - error(E_MEM,"sp_resize");
- - for ( i = A->m; i < m; i++ )
- - {
- - if ( ! (A->row[i].elt = NEW_A(MINROWLEN,row_elt)) )
- - error(E_MEM,"sp_resize");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
- - }
- - A->row[i].len = 0; A->row[i].maxlen = MINROWLEN;
- - }
- - A->m = A->max_m = m;
- - }
- -
- - /* update number of rows */
- - A->n = n;
- -
- - /* do we need to increase the size of start_idx[] and start_row[] ? */
- - if ( n > A->max_n )
- - { /* only have to update the start_idx & start_row arrays */
- - if (mem_info_is_on())
- - {
- - mem_bytes(TYPE_SPMAT,2*A->max_n*sizeof(int),
- - 2*n*sizeof(int));
- - }
- -
- - A->start_row = RENEW(A->start_row,(unsigned)n,int);
- - A->start_idx = RENEW(A->start_idx,(unsigned)n,int);
- - if ( ! A->start_row || ! A->start_idx )
- - error(E_MEM,"sp_resize");
- - A->max_n = n; /* ...and update max_n */
- -
- - return A;
- - }
- -
- - if ( n <= A->n )
- - /* make sure that all rows are truncated just before column n */
- - for ( i = 0; i < A->m; i++ )
- - {
- - r = &(A->row[i]);
- - len = sprow_idx(r,n);
- - if ( len < 0 )
- - len = -(len+2);
- - if ( len < 0 )
- - error(E_MEM,"sp_resize");
- - r->len = len;
- - }
- -
- - return A;
- -}
- -
- -
- -/* sp_compact -- removes zeros and near-zeros from a sparse matrix */
- -SPMAT *sp_compact(A,tol)
- -SPMAT *A;
- -double tol;
- -{
- - int i, idx1, idx2;
- - SPROW *r;
- - row_elt *elt1, *elt2;
- -
- - if ( ! A )
- - error(E_NULL,"sp_compact");
- - if ( tol < 0.0 )
- - error(E_RANGE,"sp_compact");
- -
- - A->flag_col = A->flag_diag = FALSE;
- -
- - for ( i = 0; i < A->m; i++ )
- - {
- - r = &(A->row[i]);
- - elt1 = elt2 = r->elt;
- - idx1 = idx2 = 0;
- - while ( idx1 < r->len )
- - {
- - /* printf("# sp_compact: idx1 = %d, idx2 = %d\n",idx1,idx2); */
- - if ( fabs(elt1->val) <= tol )
- - { idx1++; elt1++; continue; }
- - if ( elt1 != elt2 )
- - MEM_COPY(elt1,elt2,sizeof(row_elt));
- - idx1++; elt1++;
- - idx2++; elt2++;
- - }
- - r->len = idx2;
- - }
- -
- - return A;
- -}
- -
- -/* varying number of arguments */
- -
- -#ifdef ANSI_C
- -
- -/* To allocate memory to many arguments.
- - The function should be called:
- - sp_get_vars(m,n,deg,&x,&y,&z,...,NULL);
- - where
- - int m,n,deg;
- - SPMAT *x, *y, *z,...;
- - The last argument should be NULL !
- - m x n is the dimension of matrices x,y,z,...
- - returned value is equal to the number of allocated variables
- -*/
- -
- -int sp_get_vars(int m,int n,int deg,...)
- -{
- - va_list ap;
- - int i=0;
- - SPMAT **par;
- -
- - va_start(ap, deg);
- - while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
- - *par = sp_get(m,n,deg);
- - i++;
- - }
- -
- - va_end(ap);
- - return i;
- -}
- -
- -
- -/* To resize memory for many arguments.
- - The function should be called:
- - sp_resize_vars(m,n,&x,&y,&z,...,NULL);
- - where
- - int m,n;
- - SPMAT *x, *y, *z,...;
- - The last argument should be NULL !
- - m X n is the resized dimension of matrices x,y,z,...
- - returned value is equal to the number of allocated variables.
- - If one of x,y,z,.. arguments is NULL then memory is allocated to this
- - argument.
- -*/
- -
- -int sp_resize_vars(int m,int n,...)
- -{
- - va_list ap;
- - int i=0;
- - SPMAT **par;
- -
- - va_start(ap, n);
- - while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
- - *par = sp_resize(*par,m,n);
- - i++;
- - }
- -
- - va_end(ap);
- - return i;
- -}
- -
- -/* To deallocate memory for many arguments.
- - The function should be called:
- - sp_free_vars(&x,&y,&z,...,NULL);
- - where
- - SPMAT *x, *y, *z,...;
- - The last argument should be NULL !
- - There must be at least one not NULL argument.
- - returned value is equal to the number of allocated variables.
- - Returned value of x,y,z,.. is VNULL.
- -*/
- -
- -int sp_free_vars(SPMAT **va,...)
- -{
- - va_list ap;
- - int i=1;
- - SPMAT **par;
- -
- - sp_free(*va);
- - *va = (SPMAT *) NULL;
- - va_start(ap, va);
- - while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
- - sp_free(*par);
- - *par = (SPMAT *)NULL;
- - i++;
- - }
- -
- - va_end(ap);
- - return i;
- -}
- -
- -
- -#elif VARARGS
- -
- -/* To allocate memory to many arguments.
- - The function should be called:
- - sp_get_vars(m,n,deg,&x,&y,&z,...,NULL);
- - where
- - int m,n,deg;
- - SPMAT *x, *y, *z,...;
- - The last argument should be NULL !
- - m x n is the dimension of matrices x,y,z,...
- - returned value is equal to the number of allocated variables
- -*/
- -
- -int sp_get_vars(va_alist) va_dcl
- -{
- - va_list ap;
- - int i=0, m, n, deg;
- - SPMAT **par;
- -
- - va_start(ap);
- - m = va_arg(ap,int);
- - n = va_arg(ap,int);
- - deg = va_arg(ap,int);
- - while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
- - *par = sp_get(m,n,deg);
- - i++;
- - }
- -
- - va_end(ap);
- - return i;
- -}
- -
- -
- -/* To resize memory for many arguments.
- - The function should be called:
- - sp_resize_vars(m,n,&x,&y,&z,...,NULL);
- - where
- - int m,n;
- - SPMAT *x, *y, *z,...;
- - The last argument should be NULL !
- - m X n is the resized dimension of matrices x,y,z,...
- - returned value is equal to the number of allocated variables.
- - If one of x,y,z,.. arguments is NULL then memory is allocated to this
- - argument.
- -*/
- -
- -int sp_resize_vars(va_alist) va_dcl
- -{
- - va_list ap;
- - int i=0, m, n;
- - SPMAT **par;
- -
- - va_start(ap);
- - m = va_arg(ap,int);
- - n = va_arg(ap,int);
- - while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
- - *par = sp_resize(*par,m,n);
- - i++;
- - }
- -
- - va_end(ap);
- - return i;
- -}
- -
- -
- -
- -/* To deallocate memory for many arguments.
- - The function should be called:
- - sp_free_vars(&x,&y,&z,...,NULL);
- - where
- - SPMAT *x, *y, *z,...;
- - The last argument should be NULL !
- - There must be at least one not NULL argument.
- - returned value is equal to the number of allocated variables.
- - Returned value of x,y,z,.. is VNULL.
- -*/
- -
- -int sp_free_vars(va_alist) va_dcl
- -{
- - va_list ap;
- - int i=0;
- - SPMAT **par;
- -
- - va_start(ap);
- - while (par = va_arg(ap,SPMAT **)) { /* NULL ends the list*/
- - sp_free(*par);
- - *par = (SPMAT *)NULL;
- - i++;
- - }
- -
- - va_end(ap);
- - return i;
- -}
- -
- -
- -
- -#endif
- -
- //GO.SYSIN DD sparse.c
- echo sprow.c 1>&2
- sed >sprow.c <<'//GO.SYSIN DD sprow.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -/*
- - Sparse rows package
- - See also: sparse.h, matrix.h
- - */
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include <stdlib.h>
- -#include "sparse.h"
- -
- -
- -static char rcsid[] = "$Id: sprow.c,v 1.1 1994/01/13 05:35:36 des Exp $";
- -
- -#define MINROWLEN 10
- -
- -
- -/* sprow_dump - prints relevant information about the sparse row r */
- -
- -void sprow_dump(fp,r)
- -FILE *fp;
- -SPROW *r;
- -{
- - int j_idx;
- - row_elt *elts;
- -
- - fprintf(fp,"SparseRow dump:\n");
- - if ( ! r )
- - { fprintf(fp,"*** NULL row ***\n"); return; }
- -
- - fprintf(fp,"row: len = %d, maxlen = %d, diag idx = %d\n",
- - r->len,r->maxlen,r->diag);
- - fprintf(fp,"element list @ 0x%lx\n",(long)(r->elt));
- - if ( ! r->elt )
- - {
- - fprintf(fp,"*** NULL element list ***\n");
- - return;
- - }
- - elts = r->elt;
- - for ( j_idx = 0; j_idx < r->len; j_idx++, elts++ )
- - fprintf(fp,"Col: %d, Val: %g, nxt_row = %d, nxt_idx = %d\n",
- - elts->col,elts->val,elts->nxt_row,elts->nxt_idx);
- - fprintf(fp,"\n");
- -}
- -
- -
- -/* sprow_idx -- get index into row for a given column in a given row
- - -- return -1 on error
- - -- return -(idx+2) where idx is index to insertion point */
- -int sprow_idx(r,col)
- -SPROW *r;
- -int col;
- -{
- - register int lo, hi, mid;
- - int tmp;
- - register row_elt *r_elt;
- -
- - /*******************************************
- - if ( r == (SPROW *)NULL )
- - return -1;
- - if ( col < 0 )
- - return -1;
- - *******************************************/
- -
- - r_elt = r->elt;
- - if ( r->len <= 0 )
- - return -2;
- -
- - /* try the hint */
- - /* if ( hint >= 0 && hint < r->len && r_elt[hint].col == col )
- - return hint; */
- -
- - /* otherwise use binary search... */
- - /* code from K&R Ch. 6, p. 125 */
- - lo = 0; hi = r->len - 1; mid = lo;
- - while ( lo <= hi )
- - {
- - mid = (hi + lo)/2;
- - if ( (tmp=r_elt[mid].col-col) > 0 )
- - hi = mid-1;
- - else if ( tmp < 0 )
- - lo = mid+1;
- - else /* tmp == 0 */
- - return mid;
- - }
- - tmp = r_elt[mid].col - col;
- -
- - if ( tmp > 0 )
- - return -(mid+2); /* insert at mid */
- - else /* tmp < 0 */
- - return -(mid+3); /* insert at mid+1 */
- -}
- -
- -
- -/* sprow_get -- gets, initialises and returns a SPROW structure
- - -- max. length is maxlen */
- -SPROW *sprow_get(maxlen)
- -int maxlen;
- -{
- - SPROW *r;
- -
- - if ( maxlen < 0 )
- - error(E_NEG,"sprow_get");
- -
- - r = NEW(SPROW);
- - if ( ! r )
- - error(E_MEM,"sprow_get");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPROW,0,sizeof(SPROW));
- - mem_numvar(TYPE_SPROW,1);
- - }
- - r->elt = NEW_A(maxlen,row_elt);
- - if ( ! r->elt )
- - error(E_MEM,"sprow_get");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPROW,0,maxlen*sizeof(row_elt));
- - }
- - r->len = 0;
- - r->maxlen = maxlen;
- - r->diag = -1;
- -
- - return r;
- -}
- -
- -
- -/* sprow_xpd -- expand row by means of realloc()
- - -- type must be TYPE_SPMAT if r is a row of a SPMAT structure,
- - otherwise it must be TYPE_SPROW
- - -- returns r */
- -SPROW *sprow_xpd(r,n,type)
- -SPROW *r;
- -int n,type;
- -{
- - int newlen;
- -
- - if ( ! r ) {
- - r = NEW(SPROW);
- - if (! r )
- - error(E_MEM,"sprow_xpd");
- - else if ( mem_info_is_on()) {
- - if (type != TYPE_SPMAT && type != TYPE_SPROW)
- - warning(WARN_WRONG_TYPE,"sprow_xpd");
- - mem_bytes(type,0,sizeof(SPROW));
- - if (type == TYPE_SPROW)
- - mem_numvar(type,1);
- - }
- - }
- -
- - if ( ! r->elt )
- - {
- - r->elt = NEW_A((unsigned)n,row_elt);
- - if ( ! r->elt )
- - error(E_MEM,"sprow_xpd");
- - else if (mem_info_is_on()) {
- - mem_bytes(type,0,n*sizeof(row_elt));
- - }
- - r->len = 0;
- - r->maxlen = n;
- - return r;
- - }
- - if ( n <= r->len )
- - newlen = max(2*r->len + 1,MINROWLEN);
- - else
- - newlen = n;
- - if ( newlen <= r->maxlen )
- - {
- - MEM_ZERO((char *)(&(r->elt[r->len])),
- - (newlen-r->len)*sizeof(row_elt));
- - r->len = newlen;
- - }
- - else
- - {
- - if (mem_info_is_on()) {
- - mem_bytes(type,r->maxlen*sizeof(row_elt),
- - newlen*sizeof(row_elt));
- - }
- - r->elt = RENEW(r->elt,newlen,row_elt);
- - if ( ! r->elt )
- - error(E_MEM,"sprow_xpd");
- - r->maxlen = newlen;
- - r->len = newlen;
- - }
- -
- - return r;
- -}
- -
- -/* sprow_resize -- resize a SPROW variable by means of realloc()
- - -- n is a new size
- - -- returns r */
- -SPROW *sprow_resize(r,n,type)
- -SPROW *r;
- -int n,type;
- -{
- - if (n < 0)
- - error(E_NEG,"sprow_resize");
- -
- - if ( ! r )
- - return sprow_get(n);
- -
- - if (n == r->len)
- - return r;
- -
- - if ( ! r->elt )
- - {
- - r->elt = NEW_A((unsigned)n,row_elt);
- - if ( ! r->elt )
- - error(E_MEM,"sprow_resize");
- - else if (mem_info_is_on()) {
- - mem_bytes(type,0,n*sizeof(row_elt));
- - }
- - r->maxlen = r->len = n;
- - return r;
- - }
- -
- - if ( n <= r->maxlen )
- - r->len = n;
- - else
- - {
- - if (mem_info_is_on()) {
- - mem_bytes(type,r->maxlen*sizeof(row_elt),
- - n*sizeof(row_elt));
- - }
- - r->elt = RENEW(r->elt,n,row_elt);
- - if ( ! r->elt )
- - error(E_MEM,"sprow_resize");
- - r->maxlen = r->len = n;
- - }
- -
- - return r;
- -}
- -
- -
- -/* release a row of a matrix */
- -int sprow_free(r)
- -SPROW *r;
- -{
- - if ( ! r )
- - return -1;
- -
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPROW,sizeof(SPROW),0);
- - mem_numvar(TYPE_SPROW,-1);
- - }
- -
- - if ( r->elt )
- - {
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPROW,r->maxlen*sizeof(row_elt),0);
- - }
- - free((char *)r->elt);
- - }
- - free((char *)r);
- - return 0;
- -}
- -
- -
- -/* sprow_merge -- merges r1 and r2 into r_out
- - -- cannot be done in-situ
- - -- type must be SPMAT or SPROW depending on
- - whether r_out is a row of a SPMAT structure
- - or a SPROW variable
- - -- returns r_out */
- -SPROW *sprow_merge(r1,r2,r_out,type)
- -SPROW *r1, *r2, *r_out;
- -int type;
- -{
- - int idx1, idx2, idx_out, len1, len2, len_out;
- - row_elt *elt1, *elt2, *elt_out;
- -
- - if ( ! r1 || ! r2 )
- - error(E_NULL,"sprow_merge");
- - if ( ! r_out )
- - r_out = sprow_get(MINROWLEN);
- - if ( r1 == r_out || r2 == r_out )
- - error(E_INSITU,"sprow_merge");
- -
- - /* Initialise */
- - len1 = r1->len; len2 = r2->len; len_out = r_out->maxlen;
- - idx1 = idx2 = idx_out = 0;
- - elt1 = r1->elt; elt2 = r2->elt; elt_out = r_out->elt;
- -
- - while ( idx1 < len1 || idx2 < len2 )
- - {
- - if ( idx_out >= len_out )
- - { /* r_out is too small */
- - r_out->len = idx_out;
- - r_out = sprow_xpd(r_out,0,type);
- - len_out = r_out->len;
- - elt_out = &(r_out->elt[idx_out]);
- - }
- - if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
- - {
- - elt_out->col = elt1->col;
- - elt_out->val = elt1->val;
- - if ( elt1->col == elt2->col && idx2 < len2 )
- - { elt2++; idx2++; }
- - elt1++; idx1++;
- - }
- - else
- - {
- - elt_out->col = elt2->col;
- - elt_out->val = elt2->val;
- - elt2++; idx2++;
- - }
- - elt_out++; idx_out++;
- - }
- - r_out->len = idx_out;
- -
- - return r_out;
- -}
- -
- -/* sprow_copy -- copies r1 and r2 into r_out
- - -- cannot be done in-situ
- - -- type must be SPMAT or SPROW depending on
- - whether r_out is a row of a SPMAT structure
- - or a SPROW variable
- - -- returns r_out */
- -SPROW *sprow_copy(r1,r2,r_out,type)
- -SPROW *r1, *r2, *r_out;
- -int type;
- -{
- - int idx1, idx2, idx_out, len1, len2, len_out;
- - row_elt *elt1, *elt2, *elt_out;
- -
- - if ( ! r1 || ! r2 )
- - error(E_NULL,"sprow_copy");
- - if ( ! r_out )
- - r_out = sprow_get(MINROWLEN);
- - if ( r1 == r_out || r2 == r_out )
- - error(E_INSITU,"sprow_copy");
- -
- - /* Initialise */
- - len1 = r1->len; len2 = r2->len; len_out = r_out->maxlen;
- - idx1 = idx2 = idx_out = 0;
- - elt1 = r1->elt; elt2 = r2->elt; elt_out = r_out->elt;
- -
- - while ( idx1 < len1 || idx2 < len2 )
- - {
- - while ( idx_out >= len_out )
- - { /* r_out is too small */
- - r_out->len = idx_out;
- - r_out = sprow_xpd(r_out,0,type);
- - len_out = r_out->maxlen;
- - elt_out = &(r_out->elt[idx_out]);
- - }
- - if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
- - {
- - elt_out->col = elt1->col;
- - elt_out->val = elt1->val;
- - if ( elt1->col == elt2->col && idx2 < len2 )
- - { elt2++; idx2++; }
- - elt1++; idx1++;
- - }
- - else
- - {
- - elt_out->col = elt2->col;
- - elt_out->val = 0.0;
- - elt2++; idx2++;
- - }
- - elt_out++; idx_out++;
- - }
- - r_out->len = idx_out;
- -
- - return r_out;
- -}
- -
- -/* sprow_mltadd -- sets r_out <- r1 + alpha.r2
- - -- cannot be in situ
- - -- only for columns j0, j0+1, ...
- - -- type must be SPMAT or SPROW depending on
- - whether r_out is a row of a SPMAT structure
- - or a SPROW variable
- - -- returns r_out */
- -SPROW *sprow_mltadd(r1,r2,alpha,j0,r_out,type)
- -SPROW *r1, *r2, *r_out;
- -double alpha;
- -int j0, type;
- -{
- - int idx1, idx2, idx_out, len1, len2, len_out;
- - row_elt *elt1, *elt2, *elt_out;
- -
- - if ( ! r1 || ! r2 )
- - error(E_NULL,"sprow_mltadd");
- - if ( r1 == r_out || r2 == r_out )
- - error(E_INSITU,"sprow_mltadd");
- - if ( j0 < 0 )
- - error(E_BOUNDS,"sprow_mltadd");
- - if ( ! r_out )
- - r_out = sprow_get(MINROWLEN);
- -
- - /* Initialise */
- - len1 = r1->len; len2 = r2->len; len_out = r_out->maxlen;
- - /* idx1 = idx2 = idx_out = 0; */
- - idx1 = sprow_idx(r1,j0);
- - idx2 = sprow_idx(r2,j0);
- - idx_out = sprow_idx(r_out,j0);
- - idx1 = (idx1 < 0) ? -(idx1+2) : idx1;
- - idx2 = (idx2 < 0) ? -(idx2+2) : idx2;
- - idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
- - elt1 = &(r1->elt[idx1]);
- - elt2 = &(r2->elt[idx2]);
- - elt_out = &(r_out->elt[idx_out]);
- -
- - while ( idx1 < len1 || idx2 < len2 )
- - {
- - if ( idx_out >= len_out )
- - { /* r_out is too small */
- - r_out->len = idx_out;
- - r_out = sprow_xpd(r_out,0,type);
- - len_out = r_out->maxlen;
- - elt_out = &(r_out->elt[idx_out]);
- - }
- - if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
- - {
- - elt_out->col = elt1->col;
- - elt_out->val = elt1->val;
- - if ( idx2 < len2 && elt1->col == elt2->col )
- - {
- - elt_out->val += alpha*elt2->val;
- - elt2++; idx2++;
- - }
- - elt1++; idx1++;
- - }
- - else
- - {
- - elt_out->col = elt2->col;
- - elt_out->val = alpha*elt2->val;
- - elt2++; idx2++;
- - }
- - elt_out++; idx_out++;
- - }
- - r_out->len = idx_out;
- -
- - return r_out;
- -}
- -
- -/* sprow_add -- sets r_out <- r1 + r2
- - -- cannot be in situ
- - -- only for columns j0, j0+1, ...
- - -- type must be SPMAT or SPROW depending on
- - whether r_out is a row of a SPMAT structure
- - or a SPROW variable
- - -- returns r_out */
- -SPROW *sprow_add(r1,r2,j0,r_out,type)
- -SPROW *r1, *r2, *r_out;
- -int j0, type;
- -{
- - int idx1, idx2, idx_out, len1, len2, len_out;
- - row_elt *elt1, *elt2, *elt_out;
- -
- - if ( ! r1 || ! r2 )
- - error(E_NULL,"sprow_add");
- - if ( r1 == r_out || r2 == r_out )
- - error(E_INSITU,"sprow_add");
- - if ( j0 < 0 )
- - error(E_BOUNDS,"sprow_add");
- - if ( ! r_out )
- - r_out = sprow_get(MINROWLEN);
- -
- - /* Initialise */
- - len1 = r1->len; len2 = r2->len; len_out = r_out->maxlen;
- - /* idx1 = idx2 = idx_out = 0; */
- - idx1 = sprow_idx(r1,j0);
- - idx2 = sprow_idx(r2,j0);
- - idx_out = sprow_idx(r_out,j0);
- - idx1 = (idx1 < 0) ? -(idx1+2) : idx1;
- - idx2 = (idx2 < 0) ? -(idx2+2) : idx2;
- - idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
- - elt1 = &(r1->elt[idx1]);
- - elt2 = &(r2->elt[idx2]);
- - elt_out = &(r_out->elt[idx_out]);
- -
- - while ( idx1 < len1 || idx2 < len2 )
- - {
- - if ( idx_out >= len_out )
- - { /* r_out is too small */
- - r_out->len = idx_out;
- - r_out = sprow_xpd(r_out,0,type);
- - len_out = r_out->maxlen;
- - elt_out = &(r_out->elt[idx_out]);
- - }
- - if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
- - {
- - elt_out->col = elt1->col;
- - elt_out->val = elt1->val;
- - if ( idx2 < len2 && elt1->col == elt2->col )
- - {
- - elt_out->val += elt2->val;
- - elt2++; idx2++;
- - }
- - elt1++; idx1++;
- - }
- - else
- - {
- - elt_out->col = elt2->col;
- - elt_out->val = elt2->val;
- - elt2++; idx2++;
- - }
- - elt_out++; idx_out++;
- - }
- - r_out->len = idx_out;
- -
- - return r_out;
- -}
- -
- -/* sprow_sub -- sets r_out <- r1 - r2
- - -- cannot be in situ
- - -- only for columns j0, j0+1, ...
- - -- type must be SPMAT or SPROW depending on
- - whether r_out is a row of a SPMAT structure
- - or a SPROW variable
- - -- returns r_out */
- -SPROW *sprow_sub(r1,r2,j0,r_out,type)
- -SPROW *r1, *r2, *r_out;
- -int j0, type;
- -{
- - int idx1, idx2, idx_out, len1, len2, len_out;
- - row_elt *elt1, *elt2, *elt_out;
- -
- - if ( ! r1 || ! r2 )
- - error(E_NULL,"sprow_sub");
- - if ( r1 == r_out || r2 == r_out )
- - error(E_INSITU,"sprow_sub");
- - if ( j0 < 0 )
- - error(E_BOUNDS,"sprow_sub");
- - if ( ! r_out )
- - r_out = sprow_get(MINROWLEN);
- -
- - /* Initialise */
- - len1 = r1->len; len2 = r2->len; len_out = r_out->maxlen;
- - /* idx1 = idx2 = idx_out = 0; */
- - idx1 = sprow_idx(r1,j0);
- - idx2 = sprow_idx(r2,j0);
- - idx_out = sprow_idx(r_out,j0);
- - idx1 = (idx1 < 0) ? -(idx1+2) : idx1;
- - idx2 = (idx2 < 0) ? -(idx2+2) : idx2;
- - idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
- - elt1 = &(r1->elt[idx1]);
- - elt2 = &(r2->elt[idx2]);
- - elt_out = &(r_out->elt[idx_out]);
- -
- - while ( idx1 < len1 || idx2 < len2 )
- - {
- - if ( idx_out >= len_out )
- - { /* r_out is too small */
- - r_out->len = idx_out;
- - r_out = sprow_xpd(r_out,0,type);
- - len_out = r_out->maxlen;
- - elt_out = &(r_out->elt[idx_out]);
- - }
- - if ( idx2 >= len2 || (idx1 < len1 && elt1->col <= elt2->col) )
- - {
- - elt_out->col = elt1->col;
- - elt_out->val = elt1->val;
- - if ( idx2 < len2 && elt1->col == elt2->col )
- - {
- - elt_out->val -= elt2->val;
- - elt2++; idx2++;
- - }
- - elt1++; idx1++;
- - }
- - else
- - {
- - elt_out->col = elt2->col;
- - elt_out->val = -elt2->val;
- - elt2++; idx2++;
- - }
- - elt_out++; idx_out++;
- - }
- - r_out->len = idx_out;
- -
- - return r_out;
- -}
- -
- -
- -/* sprow_smlt -- sets r_out <- alpha*r1
- - -- can be in situ
- - -- only for columns j0, j0+1, ...
- - -- returns r_out */
- -SPROW *sprow_smlt(r1,alpha,j0,r_out,type)
- -SPROW *r1, *r_out;
- -double alpha;
- -int j0, type;
- -{
- - int idx1, idx_out, len1;
- - row_elt *elt1, *elt_out;
- -
- - if ( ! r1 )
- - error(E_NULL,"sprow_smlt");
- - if ( j0 < 0 )
- - error(E_BOUNDS,"sprow_smlt");
- - if ( ! r_out )
- - r_out = sprow_get(MINROWLEN);
- -
- - /* Initialise */
- - len1 = r1->len;
- - idx1 = sprow_idx(r1,j0);
- - idx_out = sprow_idx(r_out,j0);
- - idx1 = (idx1 < 0) ? -(idx1+2) : idx1;
- - idx_out = (idx_out < 0) ? -(idx_out+2) : idx_out;
- - elt1 = &(r1->elt[idx1]);
- -
- - r_out = sprow_resize(r_out,idx_out+len1-idx1,type);
- - elt_out = &(r_out->elt[idx_out]);
- -
- - for ( ; idx1 < len1; elt1++,elt_out++,idx1++,idx_out++ )
- - {
- - elt_out->col = elt1->col;
- - elt_out->val = alpha*elt1->val;
- - }
- -
- - r_out->len = idx_out;
- -
- - return r_out;
- -}
- -
- -
- -/* sprow_foutput -- print a representation of r on stream fp */
- -void sprow_foutput(fp,r)
- -FILE *fp;
- -SPROW *r;
- -{
- - int i, len;
- - row_elt *e;
- -
- - if ( ! r )
- - {
- - fprintf(fp,"SparseRow: **** NULL ****\n");
- - return;
- - }
- - len = r->len;
- - fprintf(fp,"SparseRow: length: %d\n",len);
- - for ( i = 0, e = r->elt; i < len; i++, e++ )
- - fprintf(fp,"Column %d: %g, next row: %d, next index %d\n",
- - e->col, e->val, e->nxt_row, e->nxt_idx);
- -}
- -
- -
- -/* sprow_set_val -- sets the j-th column entry of the sparse row r
- - -- Note: destroys the usual column & row access paths */
- -double sprow_set_val(r,j,val)
- -SPROW *r;
- -int j;
- -double val;
- -{
- - int idx, idx2, new_len;
- -
- - if ( ! r )
- - error(E_NULL,"sprow_set_val");
- -
- - idx = sprow_idx(r,j);
- - if ( idx >= 0 )
- - { r->elt[idx].val = val; return val; }
- - /* else */ if ( idx < -1 )
- - {
- - /* shift & insert new value */
- - idx = -(idx+2); /* this is the intended insertion index */
- - if ( r->len >= r->maxlen )
- - {
- - r->len = r->maxlen;
- - new_len = max(2*r->maxlen+1,5);
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPROW,r->maxlen*sizeof(row_elt),
- - new_len*sizeof(row_elt));
- - }
- -
- - r->elt = RENEW(r->elt,new_len,row_elt);
- - if ( ! r->elt ) /* can't allocate */
- - error(E_MEM,"sprow_set_val");
- - r->maxlen = 2*r->maxlen+1;
- - }
- - for ( idx2 = r->len-1; idx2 >= idx; idx2-- )
- - MEM_COPY((char *)(&(r->elt[idx2])),
- - (char *)(&(r->elt[idx2+1])),sizeof(row_elt));
- - /************************************************************
- - if ( idx < r->len )
- - MEM_COPY((char *)(&(r->elt[idx])),(char *)(&(r->elt[idx+1])),
- - (r->len-idx)*sizeof(row_elt));
- - ************************************************************/
- - r->len++;
- - r->elt[idx].col = j;
- - r->elt[idx].nxt_row = -1;
- - r->elt[idx].nxt_idx = -1;
- - return r->elt[idx].val = val;
- - }
- - /* else -- idx == -1, error in index/matrix! */
- - return 0.0;
- -}
- -
- -
- //GO.SYSIN DD sprow.c
- echo sparseio.c 1>&2
- sed >sparseio.c <<'//GO.SYSIN DD sparseio.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -
- -/*
- - This file has the routines for sparse matrix input/output
- - It works in conjunction with sparse.c, sparse.h etc
- -*/
- -
- -#include <stdio.h>
- -#include "sparse.h"
- -
- -static char rcsid[] = "$Id: sparseio.c,v 1.4 1994/01/13 05:34:25 des Exp $";
- -
- -
- -
- -/* local variables */
- -static char line[MAXLINE];
- -
- -/* sp_foutput -- output sparse matrix A to file/stream fp */
- -void sp_foutput(fp,A)
- -FILE *fp;
- -SPMAT *A;
- -{
- - int i, j_idx, m /* , n */;
- - SPROW *rows;
- - row_elt *elts;
- -
- - fprintf(fp,"SparseMatrix: ");
- - if ( A == SMNULL )
- - {
- - fprintf(fp,"*** NULL ***\n");
- - error(E_NULL,"sp_foutput"); return;
- - }
- - fprintf(fp,"%d by %d\n",A->m,A->n);
- - m = A->m; /* n = A->n; */
- - if ( ! (rows=A->row) )
- - {
- - fprintf(fp,"*** NULL rows ***\n");
- - error(E_NULL,"sp_foutput"); return;
- - }
- -
- - for ( i = 0; i < m; i++ )
- - {
- - fprintf(fp,"row %d: ",i);
- - if ( ! (elts=rows[i].elt) )
- - {
- - fprintf(fp,"*** NULL element list ***\n");
- - continue;
- - }
- - for ( j_idx = 0; j_idx < rows[i].len; j_idx++ )
- - {
- - fprintf(fp,"%d:%-20.15g ",elts[j_idx].col,
- - elts[j_idx].val);
- - if ( j_idx % 3 == 2 && j_idx != rows[i].len-1 )
- - fprintf(fp,"\n ");
- - }
- - fprintf(fp,"\n");
- - }
- - fprintf(fp,"#\n"); /* to stop looking beyond for next entry */
- -}
- -
- -/* sp_foutput2 -- print out sparse matrix **as a dense matrix**
- - -- see output format used in matrix.h etc */
- -/******************************************************************
- -void sp_foutput2(fp,A)
- -FILE *fp;
- -SPMAT *A;
- -{
- - int cnt, i, j, j_idx;
- - SPROW *r;
- - row_elt *elt;
- -
- - if ( A == SMNULL )
- - {
- - fprintf(fp,"Matrix: *** NULL ***\n");
- - return;
- - }
- - fprintf(fp,"Matrix: %d by %d\n",A->m,A->n);
- - for ( i = 0; i < A->m; i++ )
- - {
- - fprintf(fp,"row %d:",i);
- - r = &(A->row[i]);
- - elt = r->elt;
- - cnt = j = j_idx = 0;
- - while ( j_idx < r->len || j < A->n )
- - {
- - if ( j_idx >= r->len )
- - fprintf(fp,"%14.9g ",0.0);
- - else if ( j < elt[j_idx].col )
- - fprintf(fp,"%14.9g ",0.0);
- - else
- - fprintf(fp,"%14.9g ",elt[j_idx++].val);
- - if ( cnt++ % 4 == 3 )
- - fprintf(fp,"\n");
- - j++;
- - }
- - fprintf(fp,"\n");
- - }
- -}
- -******************************************************************/
- -
- -/* sp_dump -- prints ALL relevant information about the sparse matrix A */
- -void sp_dump(fp,A)
- -FILE *fp;
- -SPMAT *A;
- -{
- - int i, j, j_idx;
- - SPROW *rows;
- - row_elt *elts;
- -
- - fprintf(fp,"SparseMatrix dump:\n");
- - if ( ! A )
- - { fprintf(fp,"*** NULL ***\n"); return; }
- - fprintf(fp,"Matrix at 0x%lx\n",(long)A);
- - fprintf(fp,"Dimensions: %d by %d\n",A->m,A->n);
- - fprintf(fp,"MaxDimensions: %d by %d\n",A->max_m,A->max_n);
- - fprintf(fp,"flag_col = %d, flag_diag = %d\n",A->flag_col,A->flag_diag);
- - fprintf(fp,"start_row @ 0x%lx:\n",(long)(A->start_row));
- - for ( j = 0; j < A->n; j++ )
- - {
- - fprintf(fp,"%d ",A->start_row[j]);
- - if ( j % 10 == 9 )
- - fprintf(fp,"\n");
- - }
- - fprintf(fp,"\n");
- - fprintf(fp,"start_idx @ 0x%lx:\n",(long)(A->start_idx));
- - for ( j = 0; j < A->n; j++ )
- - {
- - fprintf(fp,"%d ",A->start_idx[j]);
- - if ( j % 10 == 9 )
- - fprintf(fp,"\n");
- - }
- - fprintf(fp,"\n");
- - fprintf(fp,"Rows @ 0x%lx:\n",(long)(A->row));
- - if ( ! A->row )
- - { fprintf(fp,"*** NULL row ***\n"); return; }
- - rows = A->row;
- - for ( i = 0; i < A->m; i++ )
- - {
- - fprintf(fp,"row %d: len = %d, maxlen = %d, diag idx = %d\n",
- - i,rows[i].len,rows[i].maxlen,rows[i].diag);
- - fprintf(fp,"element list @ 0x%lx\n",(long)(rows[i].elt));
- - if ( ! rows[i].elt )
- - {
- - fprintf(fp,"*** NULL element list ***\n");
- - continue;
- - }
- - elts = rows[i].elt;
- - for ( j_idx = 0; j_idx < rows[i].len; j_idx++, elts++ )
- - fprintf(fp,"Col: %d, Val: %g, nxt_row = %d, nxt_idx = %d\n",
- - elts->col,elts->val,elts->nxt_row,elts->nxt_idx);
- - fprintf(fp,"\n");
- - }
- -}
- -
- -#define MAXSCRATCH 100
- -
- -/* sp_finput -- input sparse matrix from stream/file fp
- - -- uses friendly input routine if fp is a tty
- - -- uses format identical to output format otherwise */
- -SPMAT *sp_finput(fp)
- -FILE *fp;
- -{
- - int i, len, ret_val;
- - int col, curr_col, m, n, tmp, tty;
- - Real val;
- - SPMAT *A;
- - SPROW *rows;
- -
- - row_elt scratch[MAXSCRATCH];
- - /* cannot handle >= MAXSCRATCH elements in a row */
- -
- - for ( i = 0; i < MAXSCRATCH; i++ )
- - scratch[i].nxt_row = scratch[i].nxt_idx = -1;
- -
- - tty = isatty(fileno(fp));
- -
- - if ( tty )
- - {
- - fprintf(stderr,"SparseMatrix: ");
- - do {
- - fprintf(stderr,"input rows cols: ");
- - if ( ! fgets(line,MAXLINE,fp) )
- - error(E_INPUT,"sp_finput");
- - } while ( sscanf(line,"%u %u",&m,&n) != 2 );
- - A = sp_get(m,n,5);
- - rows = A->row;
- -
- - for ( i = 0; i < m; i++ )
- - {
- - fprintf(stderr,"Row %d:\n",i);
- - fprintf(stderr,"Enter <col> <val> or 'e' to end row\n");
- - curr_col = -1;
- - for ( len = 0; len < MAXSCRATCH; len++ )
- - {
- - do {
- - fprintf(stderr,"Entry %d: ",len);
- - if ( ! fgets(line,MAXLINE,fp) )
- - error(E_INPUT,"sp_finput");
- - if ( *line == 'e' || *line == 'E' )
- - break;
- -#if REAL == DOUBLE
- - } while ( sscanf(line,"%u %lf",&col,&val) != 2 ||
- -#elif REAL == FLOAT
- - } while ( sscanf(line,"%u %f",&col,&val) != 2 ||
- -#endif
- - col >= n || col <= curr_col );
- -
- - if ( *line == 'e' || *line == 'E' )
- - break;
- -
- - scratch[len].col = col;
- - scratch[len].val = val;
- - curr_col = col;
- - }
- -
- - /* Note: len = # elements in row */
- - if ( len > 5 )
- - {
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_SPMAT,
- - A->row[i].maxlen*sizeof(row_elt),
- - len*sizeof(row_elt));
- - }
- -
- - rows[i].elt = (row_elt *)realloc((char *)rows[i].elt,
- - len*sizeof(row_elt));
- - rows[i].maxlen = len;
- - }
- - MEM_COPY(scratch,rows[i].elt,len*sizeof(row_elt));
- - rows[i].len = len;
- - rows[i].diag = sprow_idx(&(rows[i]),i);
- - }
- - }
- - else /* not tty */
- - {
- - ret_val = 0;
- - skipjunk(fp);
- - fscanf(fp,"SparseMatrix:");
- - skipjunk(fp);
- - if ( (ret_val=fscanf(fp,"%u by %u",&m,&n)) != 2 )
- - error((ret_val == EOF) ? E_EOF : E_FORMAT,"sp_finput");
- - A = sp_get(m,n,5);
- -
- - /* initialise start_row */
- - for ( i = 0; i < A->n; i++ )
- - A->start_row[i] = -1;
- -
- - rows = A->row;
- - for ( i = 0; i < m; i++ )
- - {
- - /* printf("Reading row # %d\n",i); */
- - rows[i].diag = -1;
- - skipjunk(fp);
- - if ( (ret_val=fscanf(fp,"row %d :",&tmp)) != 1 ||
- - tmp != i )
- - error((ret_val == EOF) ? E_EOF : E_FORMAT,
- - "sp_finput");
- - curr_col = -1;
- - for ( len = 0; len < MAXSCRATCH; len++ )
- - {
- -#if REAL == DOUBLE
- - if ( (ret_val=fscanf(fp,"%u : %lf",&col,&val)) != 2 )
- -#elif REAL == FLOAT
- - if ( (ret_val=fscanf(fp,"%u : %f",&col,&val)) != 2 )
- -#endif
- - break;
- - if ( col <= curr_col || col >= n )
- - error(E_FORMAT,"sp_finput");
- - scratch[len].col = col;
- - scratch[len].val = val;
- - }
- - if ( ret_val == EOF )
- - error(E_EOF,"sp_finput");
- -
- - if ( len > rows[i].maxlen )
- - {
- - rows[i].elt = (row_elt *)realloc((char *)rows[i].elt,
- - len*sizeof(row_elt));
- - rows[i].maxlen = len;
- - }
- - MEM_COPY(scratch,rows[i].elt,len*sizeof(row_elt));
- - rows[i].len = len;
- - /* printf("Have read row # %d\n",i); */
- - rows[i].diag = sprow_idx(&(rows[i]),i);
- - /* printf("Have set diag index for row # %d\n",i); */
- - }
- - }
- -
- - return A;
- -}
- -
- //GO.SYSIN DD sparseio.c
- echo spchfctr.c 1>&2
- sed >spchfctr.c <<'//GO.SYSIN DD spchfctr.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -
- -/*
- - Sparse Cholesky factorisation code
- - To be used with sparse.h, sparse.c etc
- -
- -*/
- -
- -static char rcsid[] = "$Id: spchfctr.c,v 1.4 1994/01/13 05:31:32 des Exp $";
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include "matrix.h"
- -#include "sparse.h"
- -#include "sparse2.h"
- -
- -
- -#ifndef MALLOCDECL
- -#ifndef ANSI_C
- -extern char *calloc(), *realloc();
- -#endif
- -#endif
- -
- -
- -
- -/* sprow_ip -- finds the (partial) inner product of a pair of sparse rows
- - -- uses a "merging" approach & assumes column ordered rows
- - -- row indices for inner product are all < lim */
- -double sprow_ip(row1, row2, lim)
- -SPROW *row1, *row2;
- -int lim;
- -{
- - int idx1, idx2, len1, len2, tmp;
- - int sprow_idx();
- - register row_elt *elts1, *elts2;
- - register Real sum;
- -
- - elts1 = row1->elt; elts2 = row2->elt;
- - len1 = row1->len; len2 = row2->len;
- -
- - sum = 0.0;
- -
- - if ( len1 <= 0 || len2 <= 0 )
- - return 0.0;
- - if ( elts1->col >= lim || elts2->col >= lim )
- - return 0.0;
- -
- - /* use sprow_idx() to speed up inner product where one row is
- - much longer than the other */
- - idx1 = idx2 = 0;
- - if ( len1 > 2*len2 )
- - {
- - idx1 = sprow_idx(row1,elts2->col);
- - idx1 = (idx1 < 0) ? -(idx1+2) : idx1;
- - if ( idx1 < 0 )
- - error(E_UNKNOWN,"sprow_ip");
- - len1 -= idx1;
- - }
- - else if ( len2 > 2*len1 )
- - {
- - idx2 = sprow_idx(row2,elts1->col);
- - idx2 = (idx2 < 0) ? -(idx2+2) : idx2;
- - if ( idx2 < 0 )
- - error(E_UNKNOWN,"sprow_ip");
- - len2 -= idx2;
- - }
- - if ( len1 <= 0 || len2 <= 0 )
- - return 0.0;
- -
- - elts1 = &(elts1[idx1]); elts2 = &(elts2[idx2]);
- -
- -
- - for ( ; ; ) /* forever do... */
- - {
- - if ( (tmp=elts1->col-elts2->col) < 0 )
- - {
- - len1--; elts1++;
- - if ( ! len1 || elts1->col >= lim )
- - break;
- - }
- - else if ( tmp > 0 )
- - {
- - len2--; elts2++;
- - if ( ! len2 || elts2->col >= lim )
- - break;
- - }
- - else
- - {
- - sum += elts1->val * elts2->val;
- - len1--; elts1++;
- - len2--; elts2++;
- - if ( ! len1 || ! len2 ||
- - elts1->col >= lim || elts2->col >= lim )
- - break;
- - }
- - }
- -
- - return sum;
- -}
- -
- -/* sprow_sqr -- returns same as sprow_ip(row, row, lim) */
- -double sprow_sqr(row, lim)
- -SPROW *row;
- -int lim;
- -{
- - register row_elt *elts;
- - int idx, len;
- - register Real sum, tmp;
- -
- - sum = 0.0;
- - elts = row->elt; len = row->len;
- - for ( idx = 0; idx < len; idx++, elts++ )
- - {
- - if ( elts->col >= lim )
- - break;
- - tmp = elts->val;
- - sum += tmp*tmp;
- - }
- -
- - return sum;
- -}
- -
- -static int *scan_row = (int *)NULL, *scan_idx = (int *)NULL,
- - *col_list = (int *)NULL;
- -static int scan_len = 0;
- -
- -/* set_scan -- expand scan_row and scan_idx arrays
- - -- return new length */
- -int set_scan(new_len)
- -int new_len;
- -{
- - if ( new_len <= scan_len )
- - return scan_len;
- - if ( new_len <= scan_len+5 )
- - new_len += 5;
- -
- - if ( ! scan_row || ! scan_idx || ! col_list )
- - {
- - scan_row = (int *)calloc(new_len,sizeof(int));
- - scan_idx = (int *)calloc(new_len,sizeof(int));
- - col_list = (int *)calloc(new_len,sizeof(int));
- - }
- - else
- - {
- - scan_row = (int *)realloc((char *)scan_row,new_len*sizeof(int));
- - scan_idx = (int *)realloc((char *)scan_idx,new_len*sizeof(int));
- - col_list = (int *)realloc((char *)col_list,new_len*sizeof(int));
- - }
- -
- - if ( ! scan_row || ! scan_idx || ! col_list )
- - error(E_MEM,"set_scan");
- - return new_len;
- -}
- -
- -/* spCHfactor -- sparse Cholesky factorisation
- - -- only the lower triangular part of A (incl. diagonal) is used */
- -SPMAT *spCHfactor(A)
- -SPMAT *A;
- -{
- - register int i;
- - int idx, k, m, minim, n, num_scan, diag_idx, tmp1;
- - Real pivot, tmp2;
- - SPROW *r_piv, *r_op;
- - row_elt *elt_piv, *elt_op, *old_elt;
- -
- - if ( A == SMNULL )
- - error(E_NULL,"spCHfactor");
- - if ( A->m != A->n )
- - error(E_SQUARE,"spCHfactor");
- -
- - /* set up access paths if not already done so */
- - sp_col_access(A);
- - sp_diag_access(A);
- -
- - /* printf("spCHfactor() -- checkpoint 1\n"); */
- - m = A->m; n = A->n;
- - for ( k = 0; k < m; k++ )
- - {
- - r_piv = &(A->row[k]);
- - if ( r_piv->len > scan_len )
- - set_scan(r_piv->len);
- - elt_piv = r_piv->elt;
- - diag_idx = sprow_idx2(r_piv,k,r_piv->diag);
- - if ( diag_idx < 0 )
- - error(E_POSDEF,"spCHfactor");
- - old_elt = &(elt_piv[diag_idx]);
- - for ( i = 0; i < r_piv->len; i++ )
- - {
- - if ( elt_piv[i].col > k )
- - break;
- - col_list[i] = elt_piv[i].col;
- - scan_row[i] = elt_piv[i].nxt_row;
- - scan_idx[i] = elt_piv[i].nxt_idx;
- - }
- - /* printf("spCHfactor() -- checkpoint 2\n"); */
- - num_scan = i; /* number of actual entries in scan_row etc. */
- - /* printf("num_scan = %d\n",num_scan); */
- -
- - /* set diagonal entry of Cholesky factor */
- - tmp2 = elt_piv[diag_idx].val - sprow_sqr(r_piv,k);
- - if ( tmp2 <= 0.0 )
- - error(E_POSDEF,"spCHfactor");
- - elt_piv[diag_idx].val = pivot = sqrt(tmp2);
- -
- - /* now set the k-th column of the Cholesky factors */
- - /* printf("k = %d\n",k); */
- - for ( ; ; ) /* forever do... */
- - {
- - /* printf("spCHfactor() -- checkpoint 3\n"); */
- - /* find next row where something (non-trivial) happens
- - i.e. find min(scan_row) */
- - /* printf("scan_row: "); */
- - minim = n;
- - for ( i = 0; i < num_scan; i++ )
- - {
- - tmp1 = scan_row[i];
- - /* printf("%d ",tmp1); */
- - minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
- - }
- - /* printf("minim = %d\n",minim); */
- - /* printf("col_list: "); */
- -/**********************************************************************
- - for ( i = 0; i < num_scan; i++ )
- - printf("%d ",col_list[i]);
- - printf("\n");
- -**********************************************************************/
- -
- - if ( minim >= n )
- - break; /* nothing more to do for this column */
- - r_op = &(A->row[minim]);
- - elt_op = r_op->elt;
- -
- - /* set next entry in column k of Cholesky factors */
- - idx = sprow_idx2(r_op,k,scan_idx[num_scan-1]);
- - if ( idx < 0 )
- - { /* fill-in */
- - sp_set_val(A,minim,k,
- - -sprow_ip(r_piv,r_op,k)/pivot);
- - /* in case a realloc() has occurred... */
- - elt_op = r_op->elt;
- - /* now set up column access path again */
- - idx = sprow_idx2(r_op,k,-(idx+2));
- - tmp1 = old_elt->nxt_row;
- - old_elt->nxt_row = minim;
- - r_op->elt[idx].nxt_row = tmp1;
- - tmp1 = old_elt->nxt_idx;
- - old_elt->nxt_idx = idx;
- - r_op->elt[idx].nxt_idx = tmp1;
- - }
- - else
- - elt_op[idx].val = (elt_op[idx].val -
- - sprow_ip(r_piv,r_op,k))/pivot;
- -
- - /* printf("spCHfactor() -- checkpoint 4\n"); */
- -
- - /* remember current element in column k for column chain */
- - idx = sprow_idx2(r_op,k,idx);
- - old_elt = &(r_op->elt[idx]);
- -
- - /* update scan_row */
- - /* printf("spCHfactor() -- checkpoint 5\n"); */
- - /* printf("minim = %d\n",minim); */
- - for ( i = 0; i < num_scan; i++ )
- - {
- - if ( scan_row[i] != minim )
- - continue;
- - idx = sprow_idx2(r_op,col_list[i],scan_idx[i]);
- - if ( idx < 0 )
- - { scan_row[i] = -1; continue; }
- - scan_row[i] = elt_op[idx].nxt_row;
- - scan_idx[i] = elt_op[idx].nxt_idx;
- - /* printf("scan_row[%d] = %d\n",i,scan_row[i]); */
- - /* printf("scan_idx[%d] = %d\n",i,scan_idx[i]); */
- - }
- -
- - }
- - /* printf("spCHfactor() -- checkpoint 6\n"); */
- - /* sp_dump(stdout,A); */
- - /* printf("\n\n\n"); */
- - }
- -
- - return A;
- -}
- -
- -/* spCHsolve -- solve L.L^T.out=b where L is a sparse matrix,
- - -- out, b dense vectors
- - -- returns out; operation may be in-situ */
- -VEC *spCHsolve(L,b,out)
- -SPMAT *L;
- -VEC *b, *out;
- -{
- - int i, j_idx, n, scan_idx, scan_row;
- - SPROW *row;
- - row_elt *elt;
- - Real diag_val, sum, *out_ve;
- -
- - if ( L == SMNULL || b == VNULL )
- - error(E_NULL,"spCHsolve");
- - if ( L->m != L->n )
- - error(E_SQUARE,"spCHsolve");
- - if ( b->dim != L->m )
- - error(E_SIZES,"spCHsolve");
- -
- - if ( ! L->flag_col )
- - sp_col_access(L);
- - if ( ! L->flag_diag )
- - sp_diag_access(L);
- -
- - out = v_copy(b,out);
- - out_ve = out->ve;
- -
- - /* forward substitution: solve L.x=b for x */
- - n = L->n;
- - for ( i = 0; i < n; i++ )
- - {
- - sum = out_ve[i];
- - row = &(L->row[i]);
- - elt = row->elt;
- - for ( j_idx = 0; j_idx < row->len; j_idx++, elt++ )
- - {
- - if ( elt->col >= i )
- - break;
- - sum -= elt->val*out_ve[elt->col];
- - }
- - if ( row->diag >= 0 )
- - out_ve[i] = sum/(row->elt[row->diag].val);
- - else
- - error(E_SING,"spCHsolve");
- - }
- -
- - /* backward substitution: solve L^T.out = x for out */
- - for ( i = n-1; i >= 0; i-- )
- - {
- - sum = out_ve[i];
- - row = &(L->row[i]);
- - /* Note that row->diag >= 0 by above loop */
- - elt = &(row->elt[row->diag]);
- - diag_val = elt->val;
- -
- - /* scan down column */
- - scan_idx = elt->nxt_idx;
- - scan_row = elt->nxt_row;
- - while ( scan_row >= 0 /* && scan_idx >= 0 */ )
- - {
- - row = &(L->row[scan_row]);
- - elt = &(row->elt[scan_idx]);
- - sum -= elt->val*out_ve[scan_row];
- - scan_idx = elt->nxt_idx;
- - scan_row = elt->nxt_row;
- - }
- - out_ve[i] = sum/diag_val;
- - }
- -
- - return out;
- -}
- -
- -/* spICHfactor -- sparse Incomplete Cholesky factorisation
- - -- does a Cholesky factorisation assuming NO FILL-IN
- - -- as for spCHfactor(), only the lower triangular part of A is used */
- -SPMAT *spICHfactor(A)
- -SPMAT *A;
- -{
- - int k, m, n, nxt_row, nxt_idx, diag_idx;
- - Real pivot, tmp2;
- - SPROW *r_piv, *r_op;
- - row_elt *elt_piv, *elt_op;
- -
- - if ( A == SMNULL )
- - error(E_NULL,"spICHfactor");
- - if ( A->m != A->n )
- - error(E_SQUARE,"spICHfactor");
- -
- - /* set up access paths if not already done so */
- - if ( ! A->flag_col )
- - sp_col_access(A);
- - if ( ! A->flag_diag )
- - sp_diag_access(A);
- -
- - m = A->m; n = A->n;
- - for ( k = 0; k < m; k++ )
- - {
- - r_piv = &(A->row[k]);
- -
- - diag_idx = r_piv->diag;
- - if ( diag_idx < 0 )
- - error(E_POSDEF,"spICHfactor");
- -
- - elt_piv = r_piv->elt;
- -
- - /* set diagonal entry of Cholesky factor */
- - tmp2 = elt_piv[diag_idx].val - sprow_sqr(r_piv,k);
- - if ( tmp2 <= 0.0 )
- - error(E_POSDEF,"spICHfactor");
- - elt_piv[diag_idx].val = pivot = sqrt(tmp2);
- -
- - /* find next row where something (non-trivial) happens */
- - nxt_row = elt_piv[diag_idx].nxt_row;
- - nxt_idx = elt_piv[diag_idx].nxt_idx;
- -
- - /* now set the k-th column of the Cholesky factors */
- - while ( nxt_row >= 0 && nxt_idx >= 0 )
- - {
- - /* nxt_row and nxt_idx give next next row (& index)
- - of the entry to be modified */
- - r_op = &(A->row[nxt_row]);
- - elt_op = r_op->elt;
- - elt_op[nxt_idx].val = (elt_op[nxt_idx].val -
- - sprow_ip(r_piv,r_op,k))/pivot;
- -
- - nxt_row = elt_op[nxt_idx].nxt_row;
- - nxt_idx = elt_op[nxt_idx].nxt_idx;
- - }
- - }
- -
- - return A;
- -}
- -
- -
- -/* spCHsymb -- symbolic sparse Cholesky factorisation
- - -- does NOT do any floating point arithmetic; just sets up the structure
- - -- only the lower triangular part of A (incl. diagonal) is used */
- -SPMAT *spCHsymb(A)
- -SPMAT *A;
- -{
- - register int i;
- - int idx, k, m, minim, n, num_scan, diag_idx, tmp1;
- - SPROW *r_piv, *r_op;
- - row_elt *elt_piv, *elt_op, *old_elt;
- -
- - if ( A == SMNULL )
- - error(E_NULL,"spCHsymb");
- - if ( A->m != A->n )
- - error(E_SQUARE,"spCHsymb");
- -
- - /* set up access paths if not already done so */
- - if ( ! A->flag_col )
- - sp_col_access(A);
- - if ( ! A->flag_diag )
- - sp_diag_access(A);
- -
- - /* printf("spCHsymb() -- checkpoint 1\n"); */
- - m = A->m; n = A->n;
- - for ( k = 0; k < m; k++ )
- - {
- - r_piv = &(A->row[k]);
- - if ( r_piv->len > scan_len )
- - set_scan(r_piv->len);
- - elt_piv = r_piv->elt;
- - diag_idx = sprow_idx2(r_piv,k,r_piv->diag);
- - if ( diag_idx < 0 )
- - error(E_POSDEF,"spCHsymb");
- - old_elt = &(elt_piv[diag_idx]);
- - for ( i = 0; i < r_piv->len; i++ )
- - {
- - if ( elt_piv[i].col > k )
- - break;
- - col_list[i] = elt_piv[i].col;
- - scan_row[i] = elt_piv[i].nxt_row;
- - scan_idx[i] = elt_piv[i].nxt_idx;
- - }
- - /* printf("spCHsymb() -- checkpoint 2\n"); */
- - num_scan = i; /* number of actual entries in scan_row etc. */
- - /* printf("num_scan = %d\n",num_scan); */
- -
- - /* now set the k-th column of the Cholesky factors */
- - /* printf("k = %d\n",k); */
- - for ( ; ; ) /* forever do... */
- - {
- - /* printf("spCHsymb() -- checkpoint 3\n"); */
- - /* find next row where something (non-trivial) happens
- - i.e. find min(scan_row) */
- - minim = n;
- - for ( i = 0; i < num_scan; i++ )
- - {
- - tmp1 = scan_row[i];
- - /* printf("%d ",tmp1); */
- - minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
- - }
- -
- - if ( minim >= n )
- - break; /* nothing more to do for this column */
- - r_op = &(A->row[minim]);
- - elt_op = r_op->elt;
- -
- - /* set next entry in column k of Cholesky factors */
- - idx = sprow_idx2(r_op,k,scan_idx[num_scan-1]);
- - if ( idx < 0 )
- - { /* fill-in */
- - sp_set_val(A,minim,k,0.0);
- - /* in case a realloc() has occurred... */
- - elt_op = r_op->elt;
- - /* now set up column access path again */
- - idx = sprow_idx2(r_op,k,-(idx+2));
- - tmp1 = old_elt->nxt_row;
- - old_elt->nxt_row = minim;
- - r_op->elt[idx].nxt_row = tmp1;
- - tmp1 = old_elt->nxt_idx;
- - old_elt->nxt_idx = idx;
- - r_op->elt[idx].nxt_idx = tmp1;
- - }
- -
- - /* printf("spCHsymb() -- checkpoint 4\n"); */
- -
- - /* remember current element in column k for column chain */
- - idx = sprow_idx2(r_op,k,idx);
- - old_elt = &(r_op->elt[idx]);
- -
- - /* update scan_row */
- - /* printf("spCHsymb() -- checkpoint 5\n"); */
- - /* printf("minim = %d\n",minim); */
- - for ( i = 0; i < num_scan; i++ )
- - {
- - if ( scan_row[i] != minim )
- - continue;
- - idx = sprow_idx2(r_op,col_list[i],scan_idx[i]);
- - if ( idx < 0 )
- - { scan_row[i] = -1; continue; }
- - scan_row[i] = elt_op[idx].nxt_row;
- - scan_idx[i] = elt_op[idx].nxt_idx;
- - /* printf("scan_row[%d] = %d\n",i,scan_row[i]); */
- - /* printf("scan_idx[%d] = %d\n",i,scan_idx[i]); */
- - }
- -
- - }
- - /* printf("spCHsymb() -- checkpoint 6\n"); */
- - }
- -
- - return A;
- -}
- -
- -/* comp_AAT -- compute A.A^T where A is a given sparse matrix */
- -SPMAT *comp_AAT(A)
- -SPMAT *A;
- -{
- - SPMAT *AAT;
- - SPROW *r, *r2;
- - row_elt *elts, *elts2;
- - int i, idx, idx2, j, m, minim, n, num_scan, tmp1;
- - Real ip;
- -
- - if ( ! A )
- - error(E_NULL,"comp_AAT");
- - m = A->m; n = A->n;
- -
- - /* set up column access paths */
- - if ( ! A->flag_col )
- - sp_col_access(A);
- -
- - AAT = sp_get(m,m,10);
- -
- - for ( i = 0; i < m; i++ )
- - {
- - /* initialisation */
- - r = &(A->row[i]);
- - elts = r->elt;
- -
- - /* set up scan lists for this row */
- - if ( r->len > scan_len )
- - set_scan(r->len);
- - for ( j = 0; j < r->len; j++ )
- - {
- - col_list[j] = elts[j].col;
- - scan_row[j] = elts[j].nxt_row;
- - scan_idx[j] = elts[j].nxt_idx;
- - }
- - num_scan = r->len;
- -
- - /* scan down the rows for next non-zero not
- - associated with a diagonal entry */
- - for ( ; ; )
- - {
- - minim = m;
- - for ( idx = 0; idx < num_scan; idx++ )
- - {
- - tmp1 = scan_row[idx];
- - minim = ( tmp1 >= 0 && tmp1 < minim ) ? tmp1 : minim;
- - }
- - if ( minim >= m )
- - break;
- - r2 = &(A->row[minim]);
- - if ( minim > i )
- - {
- - ip = sprow_ip(r,r2,n);
- - sp_set_val(AAT,minim,i,ip);
- - sp_set_val(AAT,i,minim,ip);
- - }
- - /* update scan entries */
- - elts2 = r2->elt;
- - for ( idx = 0; idx < num_scan; idx++ )
- - {
- - if ( scan_row[idx] != minim || scan_idx[idx] < 0 )
- - continue;
- - idx2 = scan_idx[idx];
- - scan_row[idx] = elts2[idx2].nxt_row;
- - scan_idx[idx] = elts2[idx2].nxt_idx;
- - }
- - }
- -
- - /* set the diagonal entry */
- - sp_set_val(AAT,i,i,sprow_sqr(r,n));
- - }
- -
- - return AAT;
- -}
- -
- //GO.SYSIN DD spchfctr.c
- echo splufctr.c 1>&2
- sed >splufctr.c <<'//GO.SYSIN DD splufctr.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Stewart & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -
- -/*
- - Sparse LU factorisation
- - See also: sparse.[ch] etc for details about sparse matrices
- -*/
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include "sparse2.h"
- -
- -
- -
- -/* Macro for speedup */
- -/* #define sprow_idx2(r,c,hint) \
- - ( ( (hint) >= 0 && (r)->elt[hint].col == (c)) ? hint : sprow_idx((r),(c)) ) */
- -
- -
- -/* spLUfactor -- sparse LU factorisation with pivoting
- - -- uses partial pivoting and Markowitz criterion
- - |a[p][k]| >= alpha * max_i |a[i][k]|
- - -- creates fill-in as needed
- - -- in situ factorisation */
- -SPMAT *spLUfactor(A,px,alpha)
- -SPMAT *A;
- -PERM *px;
- -double alpha;
- -{
- - int i, best_i, k, idx, len, best_len, m, n;
- - SPROW *r, *r_piv, tmp_row;
- - static SPROW *merge = (SPROW *)NULL;
- - Real max_val, tmp;
- - static VEC *col_vals=VNULL;
- -
- - if ( ! A || ! px )
- - error(E_NULL,"spLUfctr");
- - if ( alpha <= 0.0 || alpha > 1.0 )
- - error(E_RANGE,"alpha in spLUfctr");
- - if ( px->size <= A->m )
- - px = px_resize(px,A->m);
- - px_ident(px);
- - col_vals = v_resize(col_vals,A->m);
- - MEM_STAT_REG(col_vals,TYPE_VEC);
- -
- - m = A->m; n = A->n;
- - if ( ! A->flag_col )
- - sp_col_access(A);
- - if ( ! A->flag_diag )
- - sp_diag_access(A);
- - A->flag_col = A->flag_diag = FALSE;
- - if ( ! merge ) {
- - merge = sprow_get(20);
- - MEM_STAT_REG(merge,TYPE_SPROW);
- - }
- -
- - for ( k = 0; k < n; k++ )
- - {
- - /* find pivot row/element for partial pivoting */
- -
- - /* get first row with a non-zero entry in the k-th column */
- - max_val = 0.0;
- - for ( i = k; i < m; i++ )
- - {
- - r = &(A->row[i]);
- - idx = sprow_idx(r,k);
- - if ( idx < 0 )
- - tmp = 0.0;
- - else
- - tmp = r->elt[idx].val;
- - if ( fabs(tmp) > max_val )
- - max_val = fabs(tmp);
- - col_vals->ve[i] = tmp;
- - }
- -
- - if ( max_val == 0.0 )
- - continue;
- -
- - best_len = n+1; /* only if no possibilities */
- - best_i = -1;
- - for ( i = k; i < m; i++ )
- - {
- - tmp = fabs(col_vals->ve[i]);
- - if ( tmp == 0.0 )
- - continue;
- - if ( tmp >= alpha*max_val )
- - {
- - r = &(A->row[i]);
- - idx = sprow_idx(r,k);
- - len = (r->len) - idx;
- - if ( len < best_len )
- - {
- - best_len = len;
- - best_i = i;
- - }
- - }
- - }
- -
- - /* swap row #best_i with row #k */
- - MEM_COPY(&(A->row[best_i]),&tmp_row,sizeof(SPROW));
- - MEM_COPY(&(A->row[k]),&(A->row[best_i]),sizeof(SPROW));
- - MEM_COPY(&tmp_row,&(A->row[k]),sizeof(SPROW));
- - /* swap col_vals entries */
- - tmp = col_vals->ve[best_i];
- - col_vals->ve[best_i] = col_vals->ve[k];
- - col_vals->ve[k] = tmp;
- - px_transp(px,k,best_i);
- -
- - r_piv = &(A->row[k]);
- - for ( i = k+1; i < n; i++ )
- - {
- - /* compute and set multiplier */
- - tmp = col_vals->ve[i]/col_vals->ve[k];
- - if ( tmp != 0.0 )
- - sp_set_val(A,i,k,tmp);
- - else
- - continue;
- -
- - /* perform row operations */
- - merge->len = 0;
- - r = &(A->row[i]);
- - sprow_mltadd(r,r_piv,-tmp,k+1,merge,TYPE_SPROW);
- - idx = sprow_idx(r,k+1);
- - if ( idx < 0 )
- - idx = -(idx+2);
- - /* see if r needs expanding */
- - if ( r->maxlen < idx + merge->len )
- - sprow_xpd(r,idx+merge->len,TYPE_SPMAT);
- - r->len = idx+merge->len;
- - MEM_COPY((char *)(merge->elt),(char *)&(r->elt[idx]),
- - merge->len*sizeof(row_elt));
- - }
- - }
- -
- - return A;
- -}
- -
- -/* spLUsolve -- solve A.x = b using factored matrix A from spLUfactor()
- - -- returns x
- - -- may not be in-situ */
- -VEC *spLUsolve(A,pivot,b,x)
- -SPMAT *A;
- -PERM *pivot;
- -VEC *b, *x;
- -{
- - int i, idx, len, lim;
- - Real sum, *x_ve;
- - SPROW *r;
- - row_elt *elt;
- -
- - if ( ! A || ! b )
- - error(E_NULL,"spLUsolve");
- - if ( (pivot != PNULL && A->m != pivot->size) || A->m != b->dim )
- - error(E_SIZES,"spLUsolve");
- - if ( ! x || x->dim != A->n )
- - x = v_resize(x,A->n);
- -
- - if ( pivot != PNULL )
- - x = px_vec(pivot,b,x);
- - else
- - x = v_copy(b,x);
- -
- - x_ve = x->ve;
- - lim = min(A->m,A->n);
- - for ( i = 0; i < lim; i++ )
- - {
- - sum = x_ve[i];
- - r = &(A->row[i]);
- - len = r->len;
- - elt = r->elt;
- - for ( idx = 0; idx < len && elt->col < i; idx++, elt++ )
- - sum -= elt->val*x_ve[elt->col];
- - x_ve[i] = sum;
- - }
- -
- - for ( i = lim-1; i >= 0; i-- )
- - {
- - sum = x_ve[i];
- - r = &(A->row[i]);
- - len = r->len;
- - elt = &(r->elt[len-1]);
- - for ( idx = len-1; idx >= 0 && elt->col > i; idx--, elt-- )
- - sum -= elt->val*x_ve[elt->col];
- - if ( idx < 0 || elt->col != i || elt->val == 0.0 )
- - error(E_SING,"spLUsolve");
- - x_ve[i] = sum/elt->val;
- - }
- -
- - return x;
- -}
- -
- -/* spLUTsolve -- solve A.x = b using factored matrix A from spLUfactor()
- - -- returns x
- - -- may not be in-situ */
- -VEC *spLUTsolve(A,pivot,b,x)
- -SPMAT *A;
- -PERM *pivot;
- -VEC *b, *x;
- -{
- - int i, idx, lim, rownum;
- - Real sum, *tmp_ve;
- - /* SPROW *r; */
- - row_elt *elt;
- - static VEC *tmp=VNULL;
- -
- - if ( ! A || ! b )
- - error(E_NULL,"spLUTsolve");
- - if ( (pivot != PNULL && A->m != pivot->size) || A->m != b->dim )
- - error(E_SIZES,"spLUTsolve");
- - tmp = v_copy(b,tmp);
- - MEM_STAT_REG(tmp,TYPE_VEC);
- -
- - if ( ! A->flag_col )
- - sp_col_access(A);
- - if ( ! A->flag_diag )
- - sp_diag_access(A);
- -
- - lim = min(A->m,A->n);
- - tmp_ve = tmp->ve;
- - /* solve U^T.tmp = b */
- - for ( i = 0; i < lim; i++ )
- - {
- - sum = tmp_ve[i];
- - rownum = A->start_row[i];
- - idx = A->start_idx[i];
- - if ( rownum < 0 || idx < 0 )
- - error(E_SING,"spLUTsolve");
- - while ( rownum < i && rownum >= 0 && idx >= 0 )
- - {
- - elt = &(A->row[rownum].elt[idx]);
- - sum -= elt->val*tmp_ve[rownum];
- - rownum = elt->nxt_row;
- - idx = elt->nxt_idx;
- - }
- - if ( rownum != i )
- - error(E_SING,"spLUTsolve");
- - elt = &(A->row[rownum].elt[idx]);
- - if ( elt->val == 0.0 )
- - error(E_SING,"spLUTsolve");
- - tmp_ve[i] = sum/elt->val;
- - }
- -
- - /* now solve L^T.tmp = (old) tmp */
- - for ( i = lim-1; i >= 0; i-- )
- - {
- - sum = tmp_ve[i];
- - rownum = i;
- - idx = A->row[rownum].diag;
- - if ( idx < 0 )
- - error(E_NULL,"spLUTsolve");
- - elt = &(A->row[rownum].elt[idx]);
- - rownum = elt->nxt_row;
- - idx = elt->nxt_idx;
- - while ( rownum < lim && rownum >= 0 && idx >= 0 )
- - {
- - elt = &(A->row[rownum].elt[idx]);
- - sum -= elt->val*tmp_ve[rownum];
- - rownum = elt->nxt_row;
- - idx = elt->nxt_idx;
- - }
- - tmp_ve[i] = sum;
- - }
- -
- - if ( pivot != PNULL )
- - x = pxinv_vec(pivot,tmp,x);
- - else
- - x = v_copy(tmp,x);
- -
- - return x;
- -}
- -
- -/* spILUfactor -- sparse modified incomplete LU factorisation with
- - no pivoting
- - -- all pivot entries are ensured to be >= alpha in magnitude
- - -- setting alpha = 0 gives incomplete LU factorisation
- - -- no fill-in is generated
- - -- in situ factorisation */
- -SPMAT *spILUfactor(A,alpha)
- -SPMAT *A;
- -double alpha;
- -{
- - int i, k, idx, idx_piv, m, n, old_idx, old_idx_piv;
- - SPROW *r, *r_piv;
- - Real piv_val, tmp;
- -
- - /* printf("spILUfactor: entered\n"); */
- - if ( ! A )
- - error(E_NULL,"spILUfactor");
- - if ( alpha < 0.0 )
- - error(E_RANGE,"[alpha] in spILUfactor");
- -
- - m = A->m; n = A->n;
- - sp_diag_access(A);
- - sp_col_access(A);
- -
- - for ( k = 0; k < n; k++ )
- - {
- - /* printf("spILUfactor(l.%d): checkpoint A: k = %d\n",__LINE__,k); */
- - /* printf("spILUfactor(l.%d): A =\n", __LINE__); */
- - /* sp_output(A); */
- - r_piv = &(A->row[k]);
- - idx_piv = r_piv->diag;
- - if ( idx_piv < 0 )
- - {
- - sprow_set_val(r_piv,k,alpha);
- - idx_piv = sprow_idx(r_piv,k);
- - }
- - /* printf("spILUfactor: checkpoint B\n"); */
- - if ( idx_piv < 0 )
- - error(E_BOUNDS,"spILUfactor");
- - old_idx_piv = idx_piv;
- - piv_val = r_piv->elt[idx_piv].val;
- - /* printf("spILUfactor: checkpoint C\n"); */
- - if ( fabs(piv_val) < alpha )
- - piv_val = ( piv_val < 0.0 ) ? -alpha : alpha;
- - if ( piv_val == 0.0 ) /* alpha == 0.0 too! */
- - error(E_SING,"spILUfactor");
- -
- - /* go to next row with a non-zero in this column */
- - i = r_piv->elt[idx_piv].nxt_row;
- - old_idx = idx = r_piv->elt[idx_piv].nxt_idx;
- - while ( i >= k )
- - {
- - /* printf("spILUfactor: checkpoint D: i = %d\n",i); */
- - /* perform row operations */
- - r = &(A->row[i]);
- - /* idx = sprow_idx(r,k); */
- - /* printf("spLUfactor(l.%d) i = %d, idx = %d\n",
- - __LINE__, i, idx); */
- - if ( idx < 0 )
- - {
- - idx = r->elt[old_idx].nxt_idx;
- - i = r->elt[old_idx].nxt_row;
- - continue;
- - }
- - /* printf("spILUfactor: checkpoint E\n"); */
- - /* compute and set multiplier */
- - r->elt[idx].val = tmp = r->elt[idx].val/piv_val;
- - /* printf("spILUfactor: piv_val = %g, multiplier = %g\n",
- - piv_val, tmp); */
- - /* printf("spLUfactor(l.%d) multiplier = %g\n", __LINE__, tmp); */
- - if ( tmp == 0.0 )
- - {
- - idx = r->elt[old_idx].nxt_idx;
- - i = r->elt[old_idx].nxt_row;
- - continue;
- - }
- - /* idx = sprow_idx(r,k+1); */
- - /* if ( idx < 0 )
- - idx = -(idx+2); */
- - idx_piv++; idx++; /* now look beyond the multiplier entry */
- - /* printf("spILUfactor: checkpoint F: idx = %d, idx_piv = %d\n",
- - idx, idx_piv); */
- - while ( idx_piv < r_piv->len && idx < r->len )
- - {
- - /* printf("spILUfactor: checkpoint G: idx = %d, idx_piv = %d\n",
- - idx, idx_piv); */
- - if ( r_piv->elt[idx_piv].col < r->elt[idx].col )
- - idx_piv++;
- - else if ( r_piv->elt[idx_piv].col > r->elt[idx].col )
- - idx++;
- - else /* column numbers match */
- - {
- - /* printf("spILUfactor(l.%d) subtract %g times the ",
- - __LINE__, tmp); */
- - /* printf("(%d,%d) entry to the (%d,%d) entry\n",
- - k, r_piv->elt[idx_piv].col,
- - i, r->elt[idx].col); */
- - r->elt[idx].val -= tmp*r_piv->elt[idx_piv].val;
- - idx++; idx_piv++;
- - }
- - }
- -
- - /* bump to next row with a non-zero in column k */
- - /* printf("spILUfactor(l.%d) column = %d, row[%d] =\n",
- - __LINE__, r->elt[old_idx].col, i); */
- - /* sprow_foutput(stdout,r); */
- - i = r->elt[old_idx].nxt_row;
- - old_idx = idx = r->elt[old_idx].nxt_idx;
- - /* printf("spILUfactor(l.%d) i = %d, idx = %d\n", __LINE__, i, idx); */
- - /* and restore idx_piv to index of pivot entry */
- - idx_piv = old_idx_piv;
- - }
- - }
- - /* printf("spILUfactor: exiting\n"); */
- - return A;
- -}
- //GO.SYSIN DD splufctr.c
- echo spbkp.c 1>&2
- sed >spbkp.c <<'//GO.SYSIN DD spbkp.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -
- -/*
- - Sparse matrix Bunch--Kaufman--Parlett factorisation and solve
- - Radical revision started Thu 05th Nov 1992, 09:36:12 AM
- - to use Karen George's suggestion of leaving the the row elements unordered
- - Radical revision completed Mon 07th Dec 1992, 10:59:57 AM
- -*/
- -
- -static char rcsid[] = "$Id: spbkp.c,v 1.5 1994/01/13 05:44:35 des Exp $";
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include "matrix.h"
- -#include "sparse.h"
- -#include "sparse2.h"
- -
- -
- -#ifdef MALLOCDECL
- -#include <malloc.h>
- -#endif
- -
- -#define alpha 0.6403882032022076 /* = (1+sqrt(17))/8 */
- -
- -
- -#define btos(x) ((x) ? "TRUE" : "FALSE")
- -
- -/* assume no use of sqr() uses side-effects */
- -#define sqr(x) ((x)*(x))
- -
- -/* unord_get_idx -- returns index (encoded if entry not allocated)
- - of the element of row r with column j
- - -- uses linear search */
- -int unord_get_idx(r,j)
- -SPROW *r;
- -int j;
- -{
- - int idx;
- - row_elt *e;
- -
- - if ( ! r || ! r->elt )
- - error(E_NULL,"unord_get_idx");
- - for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
- - if ( e->col == j )
- - break;
- - if ( idx >= r->len )
- - return -(r->len+2);
- - else
- - return idx;
- -}
- -
- -/* unord_get_val -- returns value of the (i,j) entry of A
- - -- same assumptions as unord_get_idx() */
- -double unord_get_val(A,i,j)
- -SPMAT *A;
- -int i, j;
- -{
- - SPROW *r;
- - int idx;
- -
- - if ( ! A )
- - error(E_NULL,"unord_get_val");
- - if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
- - error(E_BOUNDS,"unord_get_val");
- -
- - r = &(A->row[i]);
- - idx = unord_get_idx(r,j);
- - if ( idx < 0 )
- - return 0.0;
- - else
- - return r->elt[idx].val;
- -}
- -
- -
- -/* bkp_swap_elt -- swaps the (i,j) with the (k,l) entry of sparse matrix
- - -- either or both of the entries may be unallocated */
- -static SPMAT *bkp_swap_elt(A,i1,j1,idx1,i2,j2,idx2)
- -SPMAT *A;
- -int i1, j1, idx1, i2, j2, idx2;
- -{
- - int tmp_row, tmp_idx;
- - SPROW *r1, *r2;
- - row_elt *e1, *e2;
- - Real tmp;
- -
- - if ( ! A )
- - error(E_NULL,"bkp_swap_elt");
- -
- - if ( i1 < 0 || j1 < 0 || i2 < 0 || j2 < 0 ||
- - i1 >= A->m || j1 >= A->n || i2 >= A->m || j2 >= A->n )
- - {
- - error(E_BOUNDS,"bkp_swap_elt");
- - }
- -
- - if ( i1 == i2 && j1 == j2 )
- - return A;
- - if ( idx1 < 0 && idx2 < 0 ) /* neither allocated */
- - return A;
- -
- - r1 = &(A->row[i1]); r2 = &(A->row[i2]);
- - /* if ( idx1 >= r1->len || idx2 >= r2->len )
- - error(E_BOUNDS,"bkp_swap_elt"); */
- - if ( idx1 < 0 ) /* assume not allocated */
- - {
- - idx1 = r1->len;
- - if ( idx1 >= r1->maxlen )
- - { tracecatch(sprow_xpd(r1,2*r1->maxlen+1,TYPE_SPMAT),
- - "bkp_swap_elt"); }
- - r1->len = idx1+1;
- - r1->elt[idx1].col = j1;
- - r1->elt[idx1].val = 0.0;
- - /* now patch up column access path */
- - tmp_row = -1; tmp_idx = j1;
- - chase_col(A,j1,&tmp_row,&tmp_idx,i1-1);
- -
- - if ( tmp_row < 0 )
- - {
- - r1->elt[idx1].nxt_row = A->start_row[j1];
- - r1->elt[idx1].nxt_idx = A->start_idx[j1];
- - A->start_row[j1] = i1;
- - A->start_idx[j1] = idx1;
- - }
- - else
- - {
- - row_elt *tmp_e;
- -
- - tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
- - r1->elt[idx1].nxt_row = tmp_e->nxt_row;
- - r1->elt[idx1].nxt_idx = tmp_e->nxt_idx;
- - tmp_e->nxt_row = i1;
- - tmp_e->nxt_idx = idx1;
- - }
- - }
- - else if ( r1->elt[idx1].col != j1 )
- - error(E_INTERN,"bkp_swap_elt");
- - if ( idx2 < 0 )
- - {
- - idx2 = r2->len;
- - if ( idx2 >= r2->maxlen )
- - { tracecatch(sprow_xpd(r2,2*r2->maxlen+1,TYPE_SPMAT),
- - "bkp_swap_elt"); }
- -
- - r2->len = idx2+1;
- - r2->elt[idx2].col = j2;
- - r2->elt[idx2].val = 0.0;
- - /* now patch up column access path */
- - tmp_row = -1; tmp_idx = j2;
- - chase_col(A,j2,&tmp_row,&tmp_idx,i2-1);
- - if ( tmp_row < 0 )
- - {
- - r2->elt[idx2].nxt_row = A->start_row[j2];
- - r2->elt[idx2].nxt_idx = A->start_idx[j2];
- - A->start_row[j2] = i2;
- - A->start_idx[j2] = idx2;
- - }
- - else
- - {
- - row_elt *tmp_e;
- -
- - tmp_e = &(A->row[tmp_row].elt[tmp_idx]);
- - r2->elt[idx2].nxt_row = tmp_e->nxt_row;
- - r2->elt[idx2].nxt_idx = tmp_e->nxt_idx;
- - tmp_e->nxt_row = i2;
- - tmp_e->nxt_idx = idx2;
- - }
- - }
- - else if ( r2->elt[idx2].col != j2 )
- - error(E_INTERN,"bkp_swap_elt");
- -
- - e1 = &(r1->elt[idx1]); e2 = &(r2->elt[idx2]);
- -
- - tmp = e1->val;
- - e1->val = e2->val;
- - e2->val = tmp;
- -
- - return A;
- -}
- -
- -/* bkp_bump_col -- bumps row and idx to next entry in column j */
- -row_elt *bkp_bump_col(A, j, row, idx)
- -SPMAT *A;
- -int j, *row, *idx;
- -{
- - SPROW *r;
- - row_elt *e;
- -
- - if ( *row < 0 )
- - {
- - *row = A->start_row[j];
- - *idx = A->start_idx[j];
- - }
- - else
- - {
- - r = &(A->row[*row]);
- - e = &(r->elt[*idx]);
- - if ( e->col != j )
- - error(E_INTERN,"bkp_bump_col");
- - *row = e->nxt_row;
- - *idx = e->nxt_idx;
- - }
- - if ( *row < 0 )
- - return (row_elt *)NULL;
- - else
- - return &(A->row[*row].elt[*idx]);
- -}
- -
- -/* bkp_interchange -- swap rows/cols i and j (symmetric pivot)
- - -- uses just the upper triangular part */
- -SPMAT *bkp_interchange(A, i1, i2)
- -SPMAT *A;
- -int i1, i2;
- -{
- - int tmp_row, tmp_idx;
- - int row1, row2, idx1, idx2, tmp_row1, tmp_idx1, tmp_row2, tmp_idx2;
- - SPROW *r1, *r2;
- - row_elt *e1, *e2;
- - IVEC *done_list = IVNULL;
- -
- - if ( ! A )
- - error(E_NULL,"bkp_interchange");
- - if ( i1 < 0 || i1 >= A->n || i2 < 0 || i2 >= A->n )
- - error(E_BOUNDS,"bkp_interchange");
- - if ( A->m != A->n )
- - error(E_SQUARE,"bkp_interchange");
- -
- - if ( i1 == i2 )
- - return A;
- - if ( i1 > i2 )
- - { tmp_idx = i1; i1 = i2; i2 = tmp_idx; }
- -
- - done_list = iv_resize(done_list,A->n);
- - for ( tmp_idx = 0; tmp_idx < A->n; tmp_idx++ )
- - done_list->ive[tmp_idx] = FALSE;
- - row1 = -1; idx1 = i1;
- - row2 = -1; idx2 = i2;
- - e1 = bkp_bump_col(A,i1,&row1,&idx1);
- - e2 = bkp_bump_col(A,i2,&row2,&idx2);
- -
- - while ( (row1 >= 0 && row1 < i1) || (row2 >= 0 && row2 < i1) )
- - /* Note: "row2 < i1" not "row2 < i2" as we must stop before the
- - "knee bend" */
- - {
- - if ( row1 >= 0 && row1 < i1 && ( row1 < row2 || row2 < 0 ) )
- - {
- - tmp_row1 = row1; tmp_idx1 = idx1;
- - e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
- - if ( ! done_list->ive[row1] )
- - {
- - if ( row1 == row2 )
- - bkp_swap_elt(A,row1,i1,idx1,row1,i2,idx2);
- - else
- - bkp_swap_elt(A,row1,i1,idx1,row1,i2,-1);
- - done_list->ive[row1] = TRUE;
- - }
- - row1 = tmp_row1; idx1 = tmp_idx1;
- - }
- - else if ( row2 >= 0 && row2 < i1 && ( row2 < row1 || row1 < 0 ) )
- - {
- - tmp_row2 = row2; tmp_idx2 = idx2;
- - e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
- - if ( ! done_list->ive[row2] )
- - {
- - if ( row1 == row2 )
- - bkp_swap_elt(A,row2,i1,idx1,row2,i2,idx2);
- - else
- - bkp_swap_elt(A,row2,i1,-1,row2,i2,idx2);
- - done_list->ive[row2] = TRUE;
- - }
- - row2 = tmp_row2; idx2 = tmp_idx2;
- - }
- - else if ( row1 == row2 )
- - {
- - tmp_row1 = row1; tmp_idx1 = idx1;
- - e1 = bkp_bump_col(A,i1,&tmp_row1,&tmp_idx1);
- - tmp_row2 = row2; tmp_idx2 = idx2;
- - e2 = bkp_bump_col(A,i2,&tmp_row2,&tmp_idx2);
- - if ( ! done_list->ive[row1] )
- - {
- - bkp_swap_elt(A,row1,i1,idx1,row2,i2,idx2);
- - done_list->ive[row1] = TRUE;
- - }
- - row1 = tmp_row1; idx1 = tmp_idx1;
- - row2 = tmp_row2; idx2 = tmp_idx2;
- - }
- - }
- -
- - /* ensure we are **past** the first knee */
- - while ( row2 >= 0 && row2 <= i1 )
- - e2 = bkp_bump_col(A,i2,&row2,&idx2);
- -
- - /* at/after 1st "knee bend" */
- - r1 = &(A->row[i1]);
- - idx1 = 0;
- - e1 = &(r1->elt[idx1]);
- - while ( row2 >= 0 && row2 < i2 )
- - {
- - /* used for update of e2 at end of loop */
- - tmp_row = row2; tmp_idx = idx2;
- - if ( ! done_list->ive[row2] )
- - {
- - r2 = &(A->row[row2]);
- - bkp_bump_col(A,i2,&tmp_row,&tmp_idx);
- - done_list->ive[row2] = TRUE;
- - tmp_idx1 = unord_get_idx(r1,row2);
- - tracecatch(bkp_swap_elt(A,row2,i2,idx2,i1,row2,tmp_idx1),
- - "bkp_interchange");
- - }
- -
- - /* update e1 and e2 */
- - row2 = tmp_row; idx2 = tmp_idx;
- - e2 = ( row2 >= 0 ) ? &(A->row[row2].elt[idx2]) : (row_elt *)NULL;
- - }
- -
- - idx1 = 0;
- - e1 = r1->elt;
- - while ( idx1 < r1->len )
- - {
- - if ( e1->col >= i2 || e1->col <= i1 )
- - {
- - idx1++;
- - e1++;
- - continue;
- - }
- - if ( ! done_list->ive[e1->col] )
- - {
- - tmp_idx2 = unord_get_idx(&(A->row[e1->col]),i2);
- - tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,e1->col,i2,tmp_idx2),
- - "bkp_interchange");
- - done_list->ive[e1->col] = TRUE;
- - }
- - idx1++;
- - e1++;
- - }
- -
- - /* at/after 2nd "knee bend" */
- - idx1 = 0;
- - e1 = &(r1->elt[idx1]);
- - r2 = &(A->row[i2]);
- - idx2 = 0;
- - e2 = &(r2->elt[idx2]);
- - while ( idx1 < r1->len )
- - {
- - if ( e1->col <= i2 )
- - {
- - idx1++; e1++;
- - continue;
- - }
- - if ( ! done_list->ive[e1->col] )
- - {
- - tmp_idx2 = unord_get_idx(r2,e1->col);
- - tracecatch(bkp_swap_elt(A,i1,e1->col,idx1,i2,e1->col,tmp_idx2),
- - "bkp_interchange");
- - done_list->ive[e1->col] = TRUE;
- - }
- - idx1++; e1++;
- - }
- -
- - idx2 = 0; e2 = r2->elt;
- - while ( idx2 < r2->len )
- - {
- - if ( e2->col <= i2 )
- - {
- - idx2++; e2++;
- - continue;
- - }
- - if ( ! done_list->ive[e2->col] )
- - {
- - tmp_idx1 = unord_get_idx(r1,e2->col);
- - tracecatch(bkp_swap_elt(A,i2,e2->col,idx2,i1,e2->col,tmp_idx1),
- - "bkp_interchange");
- - done_list->ive[e2->col] = TRUE;
- - }
- - idx2++; e2++;
- - }
- -
- - /* now interchange the digonal entries! */
- - idx1 = unord_get_idx(&(A->row[i1]),i1);
- - idx2 = unord_get_idx(&(A->row[i2]),i2);
- - if ( idx1 >= 0 || idx2 >= 0 )
- - {
- - tracecatch(bkp_swap_elt(A,i1,i1,idx1,i2,i2,idx2),
- - "bkp_interchange");
- - }
- -
- - return A;
- -}
- -
- -
- -/* iv_min -- returns minimum of an integer vector
- - -- sets index to the position in iv if index != NULL */
- -int iv_min(iv,index)
- -IVEC *iv;
- -int *index;
- -{
- - int i, i_min, min_val, tmp;
- -
- - if ( ! iv )
- - error(E_NULL,"iv_min");
- - if ( iv->dim <= 0 )
- - error(E_SIZES,"iv_min");
- - i_min = 0;
- - min_val = iv->ive[0];
- - for ( i = 1; i < iv->dim; i++ )
- - {
- - tmp = iv->ive[i];
- - if ( tmp < min_val )
- - {
- - min_val = tmp;
- - i_min = i;
- - }
- - }
- -
- - if ( index != (int *)NULL )
- - *index = i_min;
- -
- - return min_val;
- -}
- -
- -/* max_row_col -- returns max { |A[j][k]| : k >= i, k != j, k != l } given j
- - using symmetry and only the upper triangular part of A */
- -static double max_row_col(A,i,j,l)
- -SPMAT *A;
- -int i, j, l;
- -{
- - int row_num, idx;
- - SPROW *r;
- - row_elt *e;
- - Real max_val, tmp;
- -
- - if ( ! A )
- - error(E_NULL,"max_row_col");
- - if ( i < 0 || i > A->n || j < 0 || j >= A->n )
- - error(E_BOUNDS,"max_row_col");
- -
- - max_val = 0.0;
- -
- - idx = unord_get_idx(&(A->row[i]),j);
- - if ( idx < 0 )
- - {
- - row_num = -1; idx = j;
- - e = chase_past(A,j,&row_num,&idx,i);
- - }
- - else
- - {
- - row_num = i;
- - e = &(A->row[i].elt[idx]);
- - }
- - while ( row_num >= 0 && row_num < j )
- - {
- - if ( row_num != l )
- - {
- - tmp = fabs(e->val);
- - if ( tmp > max_val )
- - max_val = tmp;
- - }
- - e = bump_col(A,j,&row_num,&idx);
- - }
- - r = &(A->row[j]);
- - for ( idx = 0, e = r->elt; idx < r->len; idx++, e++ )
- - {
- - if ( e->col > j && e->col != l )
- - {
- - tmp = fabs(e->val);
- - if ( tmp > max_val )
- - max_val = tmp;
- - }
- - }
- -
- - return max_val;
- -}
- -
- -/* nonzeros -- counts non-zeros in A */
- -static int nonzeros(A)
- -SPMAT *A;
- -{
- - int cnt, i;
- -
- - if ( ! A )
- - return 0;
- - cnt = 0;
- - for ( i = 0; i < A->m; i++ )
- - cnt += A->row[i].len;
- -
- - return cnt;
- -}
- -
- -/* chk_col_access -- for spBKPfactor()
- - -- checks that column access path is OK */
- -int chk_col_access(A)
- -SPMAT *A;
- -{
- - int cnt_nz, j, row, idx;
- - SPROW *r;
- - row_elt *e;
- -
- - if ( ! A )
- - error(E_NULL,"chk_col_access");
- -
- - /* count nonzeros as we go down columns */
- - cnt_nz = 0;
- - for ( j = 0; j < A->n; j++ )
- - {
- - row = A->start_row[j];
- - idx = A->start_idx[j];
- - while ( row >= 0 )
- - {
- - if ( row >= A->m || idx < 0 )
- - return FALSE;
- - r = &(A->row[row]);
- - if ( idx >= r->len )
- - return FALSE;
- - e = &(r->elt[idx]);
- - if ( e->nxt_row >= 0 && e->nxt_row <= row )
- - return FALSE;
- - row = e->nxt_row;
- - idx = e->nxt_idx;
- - cnt_nz++;
- - }
- - }
- -
- - if ( cnt_nz != nonzeros(A) )
- - return FALSE;
- - else
- - return TRUE;
- -}
- -
- -/* col_cmp -- compare two columns -- for sorting rows using qsort() */
- -static int col_cmp(e1,e2)
- -row_elt *e1, *e2;
- -{
- - return e1->col - e2->col;
- -}
- -
- -/* spBKPfactor -- sparse Bunch-Kaufman-Parlett factorisation of A in-situ
- - -- A is factored into the form P'AP = MDM' where
- - P is a permutation matrix, M lower triangular and D is block
- - diagonal with blocks of size 1 or 2
- - -- P is stored in pivot; blocks[i]==i iff D[i][i] is a block */
- -SPMAT *spBKPfactor(A,pivot,blocks,tol)
- -SPMAT *A;
- -PERM *pivot, *blocks;
- -double tol;
- -{
- - int i, j, k, l, n, onebyone, r;
- - int idx, idx1, idx_piv;
- - int row_num;
- - int best_deg, best_j, best_l, best_cost, mark_cost, deg, deg_j,
- - deg_l, ignore_deg;
- - int list_idx, list_idx2, old_list_idx;
- - SPROW *row, *r_piv, *r1_piv;
- - row_elt *e, *e1;
- - Real aii, aip1, aip1i;
- - Real det, max_j, max_l, s, t;
- - static IVEC *scan_row = IVNULL, *scan_idx = IVNULL, *col_list = IVNULL,
- - *tmp_iv = IVNULL;
- - static IVEC *deg_list = IVNULL;
- - static IVEC *orig_idx = IVNULL, *orig1_idx = IVNULL;
- - static PERM *order = PNULL;
- -
- - if ( ! A || ! pivot || ! blocks )
- - error(E_NULL,"spBKPfactor");
- - if ( A->m != A->n )
- - error(E_SQUARE,"spBKPfactor");
- - if ( A->m != pivot->size || pivot->size != blocks->size )
- - error(E_SIZES,"spBKPfactor");
- - if ( tol <= 0.0 || tol > 1.0 )
- - error(E_RANGE,"spBKPfactor");
- -
- - n = A->n;
- -
- - px_ident(pivot); px_ident(blocks);
- - sp_col_access(A); sp_diag_access(A);
- - ignore_deg = FALSE;
- -
- - deg_list = iv_resize(deg_list,n);
- - order = px_resize(order,n);
- - MEM_STAT_REG(deg_list,TYPE_IVEC);
- - MEM_STAT_REG(order,TYPE_PERM);
- -
- - scan_row = iv_resize(scan_row,5);
- - scan_idx = iv_resize(scan_idx,5);
- - col_list = iv_resize(col_list,5);
- - orig_idx = iv_resize(orig_idx,5);
- - orig_idx = iv_resize(orig1_idx,5);
- - orig_idx = iv_resize(tmp_iv,5);
- - MEM_STAT_REG(scan_row,TYPE_IVEC);
- - MEM_STAT_REG(scan_idx,TYPE_IVEC);
- - MEM_STAT_REG(col_list,TYPE_IVEC);
- - MEM_STAT_REG(orig_idx,TYPE_IVEC);
- - MEM_STAT_REG(orig1_idx,TYPE_IVEC);
- - MEM_STAT_REG(tmp_iv,TYPE_IVEC);
- -
- - for ( i = 0; i < n-1; i = onebyone ? i+1 : i+2 )
- - {
- - /* now we want to use a Markowitz-style selection rule for
- - determining which rows to swap and whether to use
- - 1x1 or 2x2 pivoting */
- -
- - /* get list of degrees of nodes */
- - deg_list = iv_resize(deg_list,n-i);
- - if ( ! ignore_deg )
- - for ( j = i; j < n; j++ )
- - deg_list->ive[j-i] = 0;
- - else
- - {
- - for ( j = i; j < n; j++ )
- - deg_list->ive[j-i] = 1;
- - if ( i < n )
- - deg_list->ive[0] = 0;
- - }
- - order = px_resize(order,n-i);
- - px_ident(order);
- -
- - if ( ! ignore_deg )
- - {
- - for ( j = i; j < n; j++ )
- - {
- - /* idx = sprow_idx(&(A->row[j]),j+1); */
- - /* idx = fixindex(idx); */
- - idx = 0;
- - row = &(A->row[j]);
- - e = &(row->elt[idx]);
- - /* deg_list->ive[j-i] += row->len - idx; */
- - for ( ; idx < row->len; idx++, e++ )
- - if ( e->col >= i )
- - deg_list->ive[e->col - i]++;
- - }
- - /* now deg_list[k] == degree of node k+i */
- -
- - /* now sort them into increasing order */
- - iv_sort(deg_list,order);
- - /* now deg_list[idx] == degree of node i+order[idx] */
- - }
- -
- - /* now we can chase through the nodes in order of increasing
- - degree, picking out the ones that satisfy our stability
- - criterion */
- - list_idx = 0; r = -1;
- - best_j = best_l = -1;
- - for ( deg = 0; deg <= n; deg++ )
- - {
- - Real ajj, all, ajl;
- -
- - if ( list_idx >= deg_list->dim )
- - break; /* That's all folks! */
- - old_list_idx = list_idx;
- - while ( list_idx < deg_list->dim &&
- - deg_list->ive[list_idx] <= deg )
- - {
- - j = i+order->pe[list_idx];
- - if ( j < i )
- - continue;
- - /* can we use row/col j for a 1 x 1 pivot? */
- - /* find max_j = max_{k>=i} {|A[k][j]|,|A[j][k]|} */
- - ajj = fabs(unord_get_val(A,j,j));
- - if ( ajj == 0.0 )
- - {
- - list_idx++;
- - continue; /* can't use this for 1 x 1 pivot */
- - }
- -
- - max_j = max_row_col(A,i,j,-1);
- - if ( ajj >= tol/* *alpha */ *max_j )
- - {
- - onebyone = TRUE;
- - best_j = j;
- - best_deg = deg_list->ive[list_idx];
- - break;
- - }
- - list_idx++;
- - }
- - if ( best_j >= 0 )
- - break;
- - best_cost = 2*n; /* > any possible Markowitz cost (bound) */
- - best_j = best_l = -1;
- - list_idx = old_list_idx;
- - while ( list_idx < deg_list->dim &&
- - deg_list->ive[list_idx] <= deg )
- - {
- - j = i+order->pe[list_idx];
- - ajj = fabs(unord_get_val(A,j,j));
- - for ( list_idx2 = 0; list_idx2 < list_idx; list_idx2++ )
- - {
- - deg_j = deg;
- - deg_l = deg_list->ive[list_idx2];
- - l = i+order->pe[list_idx2];
- - if ( l < i )
- - continue;
- - /* try using rows/cols (j,l) for a 2 x 2 pivot block */
- - all = fabs(unord_get_val(A,l,l));
- - ajl = ( j > l ) ? fabs(unord_get_val(A,l,j)) :
- - fabs(unord_get_val(A,j,l));
- - det = fabs(ajj*all - ajl*ajl);
- - if ( det == 0.0 )
- - continue;
- - max_j = max_row_col(A,i,j,l);
- - max_l = max_row_col(A,i,l,j);
- - if ( tol*(all*max_j+ajl*max_l) < det &&
- - tol*(ajl*max_j+ajj*max_l) < det )
- - {
- - /* acceptably stable 2 x 2 pivot */
- - /* this is actually an overestimate of the
- - Markowitz cost for choosing (j,l) */
- - mark_cost = (ajj == 0.0) ?
- - ((all == 0.0) ? deg_j+deg_l : deg_j+2*deg_l) :
- - ((all == 0.0) ? 2*deg_j+deg_l :
- - 2*(deg_j+deg_l));
- - if ( mark_cost < best_cost )
- - {
- - onebyone = FALSE;
- - best_cost = mark_cost;
- - best_j = j;
- - best_l = l;
- - best_deg = deg_j;
- - }
- - }
- - }
- - list_idx++;
- - }
- - if ( best_j >= 0 )
- - break;
- - }
- -
- - if ( best_deg > (int)floor(0.8*(n-i)) )
- - ignore_deg = TRUE;
- -
- - /* now do actual interchanges */
- - if ( best_j >= 0 && onebyone )
- - {
- - bkp_interchange(A,i,best_j);
- - px_transp(pivot,i,best_j);
- - }
- - else if ( best_j >= 0 && best_l >= 0 && ! onebyone )
- - {
- - if ( best_j == i || best_j == i+1 )
- - {
- - if ( best_l == i || best_l == i+1 )
- - {
- - /* no pivoting, but must update blocks permutation */
- - px_transp(blocks,i,i+1);
- - goto dopivot;
- - }
- - bkp_interchange(A,(best_j == i) ? i+1 : i,best_l);
- - px_transp(pivot,(best_j == i) ? i+1 : i,best_l);
- - }
- - else if ( best_l == i || best_l == i+1 )
- - {
- - bkp_interchange(A,(best_l == i) ? i+1 : i,best_j);
- - px_transp(pivot,(best_l == i) ? i+1 : i,best_j);
- - }
- - else /* best_j & best_l outside i, i+1 */
- - {
- - if ( i != best_j )
- - {
- - bkp_interchange(A,i,best_j);
- - px_transp(pivot,i,best_j);
- - }
- - if ( i+1 != best_l )
- - {
- - bkp_interchange(A,i+1,best_l);
- - px_transp(pivot,i+1,best_l);
- - }
- - }
- - }
- - else /* can't pivot &/or nothing to pivot */
- - continue;
- -
- - /* update blocks permutation */
- - if ( ! onebyone )
- - px_transp(blocks,i,i+1);
- -
- - dopivot:
- - if ( onebyone )
- - {
- - int idx_j, idx_k, s_idx, s_idx2;
- - row_elt *e_ij, *e_ik;
- -
- - r_piv = &(A->row[i]);
- - idx_piv = unord_get_idx(r_piv,i);
- - /* if idx_piv < 0 then aii == 0 and no pivoting can be done;
- - -- this means that we should continue to the next iteration */
- - if ( idx_piv < 0 )
- - continue;
- - aii = r_piv->elt[idx_piv].val;
- - if ( aii == 0.0 )
- - continue;
- -
- - /* for ( j = i+1; j < n; j++ ) { ... pivot step ... } */
- - /* initialise scan_... etc for the 1 x 1 pivot */
- - scan_row = iv_resize(scan_row,r_piv->len);
- - scan_idx = iv_resize(scan_idx,r_piv->len);
- - col_list = iv_resize(col_list,r_piv->len);
- - orig_idx = iv_resize(orig_idx,r_piv->len);
- - row_num = i; s_idx = idx = 0;
- - e = &(r_piv->elt[idx]);
- - for ( idx = 0; idx < r_piv->len; idx++, e++ )
- - {
- - if ( e->col < i )
- - continue;
- - scan_row->ive[s_idx] = i;
- - scan_idx->ive[s_idx] = idx;
- - orig_idx->ive[s_idx] = idx;
- - col_list->ive[s_idx] = e->col;
- - s_idx++;
- - }
- - scan_row = iv_resize(scan_row,s_idx);
- - scan_idx = iv_resize(scan_idx,s_idx);
- - col_list = iv_resize(col_list,s_idx);
- - orig_idx = iv_resize(orig_idx,s_idx);
- -
- - order = px_resize(order,scan_row->dim);
- - px_ident(order);
- - iv_sort(col_list,order);
- -
- - tmp_iv = iv_resize(tmp_iv,scan_row->dim);
- - for ( idx = 0; idx < order->size; idx++ )
- - tmp_iv->ive[idx] = scan_idx->ive[order->pe[idx]];
- - iv_copy(tmp_iv,scan_idx);
- - for ( idx = 0; idx < order->size; idx++ )
- - tmp_iv->ive[idx] = scan_row->ive[order->pe[idx]];
- - iv_copy(tmp_iv,scan_row);
- - for ( idx = 0; idx < scan_row->dim; idx++ )
- - tmp_iv->ive[idx] = orig_idx->ive[order->pe[idx]];
- - iv_copy(tmp_iv,orig_idx);
- -
- - /* now do actual pivot */
- - /* for ( j = i+1; j < n-1; j++ ) .... */
- -
- - for ( s_idx = 0; s_idx < scan_row->dim; s_idx++ )
- - {
- - idx_j = orig_idx->ive[s_idx];
- - if ( idx_j < 0 )
- - error(E_INTERN,"spBKPfactor");
- - e_ij = &(r_piv->elt[idx_j]);
- - j = e_ij->col;
- - if ( j < i+1 )
- - continue;
- - scan_to(A,scan_row,scan_idx,col_list,j);
- -
- - /* compute multiplier */
- - t = e_ij->val / aii;
- -
- - /* for ( k = j; k < n; k++ ) { .... update A[j][k] .... } */
- - /* this is the row in which pivoting is done */
- - row = &(A->row[j]);
- - for ( s_idx2 = s_idx; s_idx2 < scan_row->dim; s_idx2++ )
- - {
- - idx_k = orig_idx->ive[s_idx2];
- - e_ik = &(r_piv->elt[idx_k]);
- - k = e_ik->col;
- - /* k >= j since col_list has been sorted */
- -
- - if ( scan_row->ive[s_idx2] == j )
- - { /* no fill-in -- can be done directly */
- - idx = scan_idx->ive[s_idx2];
- - /* idx = sprow_idx2(row,k,idx); */
- - row->elt[idx].val -= t*e_ik->val;
- - }
- - else
- - { /* fill-in -- insert entry & patch column */
- - int old_row, old_idx;
- - row_elt *old_e, *new_e;
- -
- - old_row = scan_row->ive[s_idx2];
- - old_idx = scan_idx->ive[s_idx2];
- - /* old_idx = sprow_idx2(&(A->row[old_row]),k,old_idx); */
- -
- - if ( old_idx < 0 )
- - error(E_INTERN,"spBKPfactor");
- - /* idx = sprow_idx(row,k); */
- - /* idx = fixindex(idx); */
- - idx = row->len;
- -
- - /* sprow_set_val(row,k,-t*e_ik->val); */
- - if ( row->len >= row->maxlen )
- - { tracecatch(sprow_xpd(row,2*row->maxlen+1,TYPE_SPMAT),
- - "spBKPfactor"); }
- -
- - row->len = idx+1;
- -
- - new_e = &(row->elt[idx]);
- - new_e->val = -t*e_ik->val;
- - new_e->col = k;
- -
- - old_e = &(A->row[old_row].elt[old_idx]);
- - new_e->nxt_row = old_e->nxt_row;
- - new_e->nxt_idx = old_e->nxt_idx;
- - old_e->nxt_row = j;
- - old_e->nxt_idx = idx;
- - }
- - }
- - e_ij->val = t;
- - }
- - }
- - else /* onebyone == FALSE */
- - { /* do 2 x 2 pivot */
- - int idx_k, idx1_k, s_idx, s_idx2;
- - int old_col;
- - row_elt *e_tmp;
- -
- - r_piv = &(A->row[i]);
- - idx_piv = unord_get_idx(r_piv,i);
- - aii = aip1i = 0.0;
- - e_tmp = r_piv->elt;
- - for ( idx_piv = 0; idx_piv < r_piv->len; idx_piv++, e_tmp++ )
- - if ( e_tmp->col == i )
- - aii = e_tmp->val;
- - else if ( e_tmp->col == i+1 )
- - aip1i = e_tmp->val;
- -
- - r1_piv = &(A->row[i+1]);
- - e_tmp = r1_piv->elt;
- - aip1 = unord_get_val(A,i+1,i+1);
- - det = aii*aip1 - aip1i*aip1i; /* Must have det < 0 */
- - if ( aii == 0.0 && aip1i == 0.0 )
- - {
- - /* error(E_RANGE,"spBKPfactor"); */
- - onebyone = TRUE;
- - continue; /* cannot pivot */
- - }
- -
- - if ( det == 0.0 )
- - {
- - if ( aii != 0.0 )
- - error(E_RANGE,"spBKPfactor");
- - onebyone = TRUE;
- - continue; /* cannot pivot */
- - }
- - aip1i = aip1i/det;
- - aii = aii/det;
- - aip1 = aip1/det;
- -
- - /* initialise scan_... etc for the 2 x 2 pivot */
- - s_idx = r_piv->len + r1_piv->len;
- - scan_row = iv_resize(scan_row,s_idx);
- - scan_idx = iv_resize(scan_idx,s_idx);
- - col_list = iv_resize(col_list,s_idx);
- - orig_idx = iv_resize(orig_idx,s_idx);
- - orig1_idx = iv_resize(orig1_idx,s_idx);
- -
- - e = r_piv->elt;
- - for ( idx = 0; idx < r_piv->len; idx++, e++ )
- - {
- - scan_row->ive[idx] = i;
- - scan_idx->ive[idx] = idx;
- - col_list->ive[idx] = e->col;
- - orig_idx->ive[idx] = idx;
- - orig1_idx->ive[idx] = -1;
- - }
- - e = r_piv->elt;
- - e1 = r1_piv->elt;
- - for ( idx = 0; idx < r1_piv->len; idx++, e1++ )
- - {
- - scan_row->ive[idx+r_piv->len] = i+1;
- - scan_idx->ive[idx+r_piv->len] = idx;
- - col_list->ive[idx+r_piv->len] = e1->col;
- - orig_idx->ive[idx+r_piv->len] = -1;
- - orig1_idx->ive[idx+r_piv->len] = idx;
- - }
- -
- - e1 = r1_piv->elt;
- - order = px_resize(order,scan_row->dim);
- - px_ident(order);
- - iv_sort(col_list,order);
- - tmp_iv = iv_resize(tmp_iv,scan_row->dim);
- - for ( idx = 0; idx < order->size; idx++ )
- - tmp_iv->ive[idx] = scan_idx->ive[order->pe[idx]];
- - iv_copy(tmp_iv,scan_idx);
- - for ( idx = 0; idx < order->size; idx++ )
- - tmp_iv->ive[idx] = scan_row->ive[order->pe[idx]];
- - iv_copy(tmp_iv,scan_row);
- - for ( idx = 0; idx < scan_row->dim; idx++ )
- - tmp_iv->ive[idx] = orig_idx->ive[order->pe[idx]];
- - iv_copy(tmp_iv,orig_idx);
- - for ( idx = 0; idx < scan_row->dim; idx++ )
- - tmp_iv->ive[idx] = orig1_idx->ive[order->pe[idx]];
- - iv_copy(tmp_iv,orig1_idx);
- -
- - s_idx = 0;
- - old_col = -1;
- - for ( idx = 0; idx < scan_row->dim; idx++ )
- - {
- - if ( col_list->ive[idx] == old_col )
- - {
- - if ( scan_row->ive[idx] == i )
- - {
- - scan_row->ive[s_idx-1] = scan_row->ive[idx];
- - scan_idx->ive[s_idx-1] = scan_idx->ive[idx];
- - col_list->ive[s_idx-1] = col_list->ive[idx];
- - orig_idx->ive[s_idx-1] = orig_idx->ive[idx];
- - orig1_idx->ive[s_idx-1] = orig1_idx->ive[idx-1];
- - }
- - else if ( idx > 0 )
- - {
- - scan_row->ive[s_idx-1] = scan_row->ive[idx-1];
- - scan_idx->ive[s_idx-1] = scan_idx->ive[idx-1];
- - col_list->ive[s_idx-1] = col_list->ive[idx-1];
- - orig_idx->ive[s_idx-1] = orig_idx->ive[idx-1];
- - orig1_idx->ive[s_idx-1] = orig1_idx->ive[idx];
- - }
- - }
- - else
- - {
- - scan_row->ive[s_idx] = scan_row->ive[idx];
- - scan_idx->ive[s_idx] = scan_idx->ive[idx];
- - col_list->ive[s_idx] = col_list->ive[idx];
- - orig_idx->ive[s_idx] = orig_idx->ive[idx];
- - orig1_idx->ive[s_idx] = orig1_idx->ive[idx];
- - s_idx++;
- - }
- - old_col = col_list->ive[idx];
- - }
- - scan_row = iv_resize(scan_row,s_idx);
- - scan_idx = iv_resize(scan_idx,s_idx);
- - col_list = iv_resize(col_list,s_idx);
- - orig_idx = iv_resize(orig_idx,s_idx);
- - orig1_idx = iv_resize(orig1_idx,s_idx);
- -
- - /* for ( j = i+2; j < n; j++ ) { .... row operation .... } */
- - for ( s_idx = 0; s_idx < scan_row->dim; s_idx++ )
- - {
- - int idx_piv, idx1_piv;
- - Real aip1j, aij, aik, aip1k;
- - row_elt *e_ik, *e_ip1k;
- -
- - j = col_list->ive[s_idx];
- - if ( j < i+2 )
- - continue;
- - tracecatch(scan_to(A,scan_row,scan_idx,col_list,j),
- - "spBKPfactor");
- -
- - idx_piv = orig_idx->ive[s_idx];
- - aij = ( idx_piv < 0 ) ? 0.0 : r_piv->elt[idx_piv].val;
- - /* aij = ( s_idx < r_piv->len ) ? r_piv->elt[s_idx].val :
- - 0.0; */
- - /* aij = sp_get_val(A,i,j); */
- - idx1_piv = orig1_idx->ive[s_idx];
- - aip1j = ( idx1_piv < 0 ) ? 0.0 : r1_piv->elt[idx1_piv].val;
- - /* aip1j = ( s_idx < r_piv->len ) ? 0.0 :
- - r1_piv->elt[s_idx-r_piv->len].val; */
- - /* aip1j = sp_get_val(A,i+1,j); */
- - s = - aip1i*aip1j + aip1*aij;
- - t = - aip1i*aij + aii*aip1j;
- -
- - /* for ( k = j; k < n; k++ ) { .... update entry .... } */
- - row = &(A->row[j]);
- - /* set idx_k and idx1_k indices */
- - s_idx2 = s_idx;
- - k = col_list->ive[s_idx2];
- - idx_k = orig_idx->ive[s_idx2];
- - idx1_k = orig1_idx->ive[s_idx2];
- -
- - while ( s_idx2 < scan_row->dim )
- - {
- - k = col_list->ive[s_idx2];
- - idx_k = orig_idx->ive[s_idx2];
- - idx1_k = orig1_idx->ive[s_idx2];
- - e_ik = ( idx_k < 0 ) ? (row_elt *)NULL :
- - &(r_piv->elt[idx_k]);
- - e_ip1k = ( idx1_k < 0 ) ? (row_elt *)NULL :
- - &(r1_piv->elt[idx1_k]);
- - aik = ( idx_k >= 0 ) ? e_ik->val : 0.0;
- - aip1k = ( idx1_k >= 0 ) ? e_ip1k->val : 0.0;
- - if ( scan_row->ive[s_idx2] == j )
- - { /* no fill-in */
- - row = &(A->row[j]);
- - /* idx = sprow_idx(row,k); */
- - idx = scan_idx->ive[s_idx2];
- - if ( idx < 0 )
- - error(E_INTERN,"spBKPfactor");
- - row->elt[idx].val -= s*aik + t*aip1k;
- - }
- - else
- - { /* fill-in -- insert entry & patch column */
- - Real tmp;
- - int old_row, old_idx;
- - row_elt *old_e, *new_e;
- -
- - tmp = - s*aik - t*aip1k;
- - if ( tmp != 0.0 )
- - {
- - row = &(A->row[j]);
- - old_row = scan_row->ive[s_idx2];
- - old_idx = scan_idx->ive[s_idx2];
- -
- - idx = row->len;
- - if ( row->len >= row->maxlen )
- - { tracecatch(sprow_xpd(row,2*row->maxlen+1,
- - TYPE_SPMAT),
- - "spBKPfactor"); }
- -
- - row->len = idx + 1;
- - /* idx = sprow_idx(row,k); */
- - new_e = &(row->elt[idx]);
- - new_e->val = tmp;
- - new_e->col = k;
- -
- - if ( old_row < 0 )
- - error(E_INTERN,"spBKPfactor");
- - /* old_idx = sprow_idx2(&(A->row[old_row]),
- - k,old_idx); */
- - old_e = &(A->row[old_row].elt[old_idx]);
- - new_e->nxt_row = old_e->nxt_row;
- - new_e->nxt_idx = old_e->nxt_idx;
- - old_e->nxt_row = j;
- - old_e->nxt_idx = idx;
- - }
- - }
- -
- - /* update idx_k, idx1_k, s_idx2 etc */
- - s_idx2++;
- - }
- -
- - /* store multipliers -- may involve fill-in (!) */
- - /* idx = sprow_idx(r_piv,j); */
- - idx = orig_idx->ive[s_idx];
- - if ( idx >= 0 )
- - {
- - r_piv->elt[idx].val = s;
- - }
- - else if ( s != 0.0 )
- - {
- - int old_row, old_idx;
- - row_elt *new_e, *old_e;
- -
- - old_row = -1; old_idx = j;
- -
- - if ( i > 0 )
- - {
- - tracecatch(chase_col(A,j,&old_row,&old_idx,i-1),
- - "spBKPfactor");
- - }
- - /* sprow_set_val(r_piv,j,s); */
- - idx = r_piv->len;
- - if ( r_piv->len >= r_piv->maxlen )
- - { tracecatch(sprow_xpd(r_piv,2*r_piv->maxlen+1,
- - TYPE_SPMAT),
- - "spBKPfactor"); }
- -
- - r_piv->len = idx + 1;
- - /* idx = sprow_idx(r_piv,j); */
- - /* if ( idx < 0 )
- - error(E_INTERN,"spBKPfactor"); */
- - new_e = &(r_piv->elt[idx]);
- - new_e->val = s;
- - new_e->col = j;
- - if ( old_row < 0 )
- - {
- - new_e->nxt_row = A->start_row[j];
- - new_e->nxt_idx = A->start_idx[j];
- - A->start_row[j] = i;
- - A->start_idx[j] = idx;
- - }
- - else
- - {
- - /* old_idx = sprow_idx2(&(A->row[old_row]),j,old_idx);*/
- - if ( old_idx < 0 )
- - error(E_INTERN,"spBKPfactor");
- - old_e = &(A->row[old_row].elt[old_idx]);
- - new_e->nxt_row = old_e->nxt_row;
- - new_e->nxt_idx = old_e->nxt_idx;
- - old_e->nxt_row = i;
- - old_e->nxt_idx = idx;
- - }
- - }
- - /* idx1 = sprow_idx(r1_piv,j); */
- - idx1 = orig1_idx->ive[s_idx];
- - if ( idx1 >= 0 )
- - {
- - r1_piv->elt[idx1].val = t;
- - }
- - else if ( t != 0.0 )
- - {
- - int old_row, old_idx;
- - row_elt *new_e, *old_e;
- -
- - old_row = -1; old_idx = j;
- - tracecatch(chase_col(A,j,&old_row,&old_idx,i),
- - "spBKPfactor");
- - /* sprow_set_val(r1_piv,j,t); */
- - idx1 = r1_piv->len;
- - if ( r1_piv->len >= r1_piv->maxlen )
- - { tracecatch(sprow_xpd(r1_piv,2*r1_piv->maxlen+1,
- - TYPE_SPMAT),
- - "spBKPfactor"); }
- -
- - r1_piv->len = idx1 + 1;
- - /* idx1 = sprow_idx(r1_piv,j); */
- - /* if ( idx < 0 )
- - error(E_INTERN,"spBKPfactor"); */
- - new_e = &(r1_piv->elt[idx1]);
- - new_e->val = t;
- - new_e->col = j;
- - if ( idx1 < 0 )
- - error(E_INTERN,"spBKPfactor");
- - new_e = &(r1_piv->elt[idx1]);
- - if ( old_row < 0 )
- - {
- - new_e->nxt_row = A->start_row[j];
- - new_e->nxt_idx = A->start_idx[j];
- - A->start_row[j] = i+1;
- - A->start_idx[j] = idx1;
- - }
- - else
- - {
- - old_idx = sprow_idx2(&(A->row[old_row]),j,old_idx);
- - if ( old_idx < 0 )
- - error(E_INTERN,"spBKPfactor");
- - old_e = &(A->row[old_row].elt[old_idx]);
- - new_e->nxt_row = old_e->nxt_row;
- - new_e->nxt_idx = old_e->nxt_idx;
- - old_e->nxt_row = i+1;
- - old_e->nxt_idx = idx1;
- - }
- - }
- - }
- - }
- - }
- -
- - /* now sort the rows arrays */
- - for ( i = 0; i < A->m; i++ )
- - qsort(A->row[i].elt,A->row[i].len,sizeof(row_elt),(int(*)())col_cmp);
- - A->flag_col = A->flag_diag = FALSE;
- -
- - return A;
- -}
- -
- -/* spBKPsolve -- solves A.x = b where A has been factored a la BKPfactor()
- - -- returns x, which is created if NULL */
- -VEC *spBKPsolve(A,pivot,block,b,x)
- -SPMAT *A;
- -PERM *pivot, *block;
- -VEC *b, *x;
- -{
- - static VEC *tmp=VNULL; /* dummy storage needed */
- - int i /* , j */, n, onebyone;
- - int row_num, idx;
- - Real a11, a12, a22, b1, b2, det, sum, *tmp_ve, tmp_diag;
- - SPROW *r;
- - row_elt *e;
- -
- - if ( ! A || ! pivot || ! block || ! b )
- - error(E_NULL,"spBKPsolve");
- - if ( A->m != A->n )
- - error(E_SQUARE,"spBKPsolve");
- - n = A->n;
- - if ( b->dim != n || pivot->size != n || block->size != n )
- - error(E_SIZES,"spBKPsolve");
- - x = v_resize(x,n);
- - tmp = v_resize(tmp,n);
- - MEM_STAT_REG(tmp,TYPE_VEC);
- -
- - tmp_ve = tmp->ve;
- -
- - if ( ! A->flag_col )
- - sp_col_access(A);
- -
- - px_vec(pivot,b,tmp);
- - /* printf("# BKPsolve: effect of pivot: tmp =\n"); v_output(tmp); */
- -
- - /* solve for lower triangular part */
- - for ( i = 0; i < n; i++ )
- - {
- - sum = tmp_ve[i];
- - if ( block->pe[i] < i )
- - {
- - /* for ( j = 0; j < i-1; j++ )
- - sum -= A_me[j][i]*tmp_ve[j]; */
- - row_num = -1; idx = i;
- - e = bump_col(A,i,&row_num,&idx);
- - while ( row_num >= 0 && row_num < i-1 )
- - {
- - sum -= e->val*tmp_ve[row_num];
- - e = bump_col(A,i,&row_num,&idx);
- - }
- - }
- - else
- - {
- - /* for ( j = 0; j < i; j++ )
- - sum -= A_me[j][i]*tmp_ve[j]; */
- - row_num = -1; idx = i;
- - e = bump_col(A,i,&row_num,&idx);
- - while ( row_num >= 0 && row_num < i )
- - {
- - sum -= e->val*tmp_ve[row_num];
- - e = bump_col(A,i,&row_num,&idx);
- - }
- - }
- - tmp_ve[i] = sum;
- - }
- -
- - /* printf("# BKPsolve: solving L part: tmp =\n"); v_output(tmp); */
- - /* solve for diagonal part */
- - for ( i = 0; i < n; i = onebyone ? i+1 : i+2 )
- - {
- - onebyone = ( block->pe[i] == i );
- - if ( onebyone )
- - {
- - /* tmp_ve[i] /= A_me[i][i]; */
- - tmp_diag = sp_get_val(A,i,i);
- - if ( tmp_diag == 0.0 )
- - error(E_SING,"spBKPsolve");
- - tmp_ve[i] /= tmp_diag;
- - }
- - else
- - {
- - a11 = sp_get_val(A,i,i);
- - a22 = sp_get_val(A,i+1,i+1);
- - a12 = sp_get_val(A,i,i+1);
- - b1 = tmp_ve[i];
- - b2 = tmp_ve[i+1];
- - det = a11*a22-a12*a12; /* < 0 : see BKPfactor() */
- - if ( det == 0.0 )
- - error(E_SING,"BKPsolve");
- - det = 1/det;
- - tmp_ve[i] = det*(a22*b1-a12*b2);
- - tmp_ve[i+1] = det*(a11*b2-a12*b1);
- - }
- - }
- -
- - /* printf("# BKPsolve: solving D part: tmp =\n"); v_output(tmp); */
- - /* solve for transpose of lower triangular part */
- - for ( i = n-2; i >= 0; i-- )
- - {
- - sum = tmp_ve[i];
- - if ( block->pe[i] > i )
- - {
- - /* onebyone is false */
- - /* for ( j = i+2; j < n; j++ )
- - sum -= A_me[i][j]*tmp_ve[j]; */
- - if ( i+2 >= n )
- - continue;
- - r = &(A->row[i]);
- - idx = sprow_idx(r,i+2);
- - idx = fixindex(idx);
- - e = &(r->elt[idx]);
- - for ( ; idx < r->len; idx++, e++ )
- - sum -= e->val*tmp_ve[e->col];
- - }
- - else /* onebyone */
- - {
- - /* for ( j = i+1; j < n; j++ )
- - sum -= A_me[i][j]*tmp_ve[j]; */
- - r = &(A->row[i]);
- - idx = sprow_idx(r,i+1);
- - idx = fixindex(idx);
- - e = &(r->elt[idx]);
- - for ( ; idx < r->len; idx++, e++ )
- - sum -= e->val*tmp_ve[e->col];
- - }
- - tmp_ve[i] = sum;
- - }
- -
- - /* printf("# BKPsolve: solving L^T part: tmp =\n");v_output(tmp); */
- - /* and do final permutation */
- - x = pxinv_vec(pivot,tmp,x);
- -
- - return x;
- -}
- -
- -
- -
- //GO.SYSIN DD spbkp.c
- echo spswap.c 1>&2
- sed >spswap.c <<'//GO.SYSIN DD spswap.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -
- -/*
- - Sparse matrix swap and permutation routines
- - Modified Mon 09th Nov 1992, 08:50:54 PM
- - to use Karen George's suggestion to use unordered rows
- -*/
- -
- -static char rcsid[] = "$Id: spswap.c,v 1.3 1994/01/13 05:44:43 des Exp $";
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include "matrix.h"
- -#include "sparse.h"
- -#include "sparse2.h"
- -
- -
- -#define btos(x) ((x) ? "TRUE" : "FALSE")
- -
- -/* scan_to -- updates scan (int) vectors to point to the last row in each
- - column with row # <= max_row, if any */
- -void scan_to(A, scan_row, scan_idx, col_list, max_row)
- -SPMAT *A;
- -IVEC *scan_row, *scan_idx, *col_list;
- -int max_row;
- -{
- - int col, idx, j_idx, row_num;
- - SPROW *r;
- - row_elt *e;
- -
- - if ( ! A || ! scan_row || ! scan_idx || ! col_list )
- - error(E_NULL,"scan_to");
- - if ( scan_row->dim != scan_idx->dim || scan_idx->dim != col_list->dim )
- - error(E_SIZES,"scan_to");
- -
- - if ( max_row < 0 )
- - return;
- -
- - if ( ! A->flag_col )
- - sp_col_access(A);
- -
- - for ( j_idx = 0; j_idx < scan_row->dim; j_idx++ )
- - {
- - row_num = scan_row->ive[j_idx];
- - idx = scan_idx->ive[j_idx];
- - col = col_list->ive[j_idx];
- -
- - if ( col < 0 || col >= A->n )
- - error(E_BOUNDS,"scan_to");
- - if ( row_num < 0 )
- - {
- - idx = col;
- - continue;
- - }
- - r = &(A->row[row_num]);
- - if ( idx < 0 )
- - error(E_INTERN,"scan_to");
- - e = &(r->elt[idx]);
- - if ( e->col != col )
- - error(E_INTERN,"scan_to");
- - if ( idx < 0 )
- - {
- - printf("scan_to: row_num = %d, idx = %d, col = %d\n",
- - row_num, idx, col);
- - error(E_INTERN,"scan_to");
- - }
- - /* if ( e->nxt_row <= max_row )
- - chase_col(A, col, &row_num, &idx, max_row); */
- - while ( e->nxt_row >= 0 && e->nxt_row <= max_row )
- - {
- - row_num = e->nxt_row;
- - idx = e->nxt_idx;
- - e = &(A->row[row_num].elt[idx]);
- - }
- -
- - /* printf("scan_to: computed j_idx = %d, row_num = %d, idx = %d\n",
- - j_idx, row_num, idx); */
- - scan_row->ive[j_idx] = row_num;
- - scan_idx->ive[j_idx] = idx;
- - }
- -}
- -
- -/* patch_col -- patches column access paths for fill-in */
- -void patch_col(A, col, old_row, old_idx, row_num, idx)
- -SPMAT *A;
- -int col, old_row, old_idx, row_num, idx;
- -{
- - SPROW *r;
- - row_elt *e;
- -
- - if ( old_row >= 0 )
- - {
- - r = &(A->row[old_row]);
- - old_idx = sprow_idx2(r,col,old_idx);
- - e = &(r->elt[old_idx]);
- - e->nxt_row = row_num;
- - e->nxt_idx = idx;
- - }
- - else
- - {
- - A->start_row[col] = row_num;
- - A->start_idx[col] = idx;
- - }
- -}
- -
- -/* chase_col -- chases column access path in column col, starting with
- - row_num and idx, to find last row # in this column <= max_row
- - -- row_num is returned; idx is also set by this routine
- - -- assumes that the column access paths (possibly without the
- - nxt_idx fields) are set up */
- -row_elt *chase_col(A, col, row_num, idx, max_row)
- -SPMAT *A;
- -int col, *row_num, *idx, max_row;
- -{
- - int old_idx, old_row, tmp_idx, tmp_row;
- - SPROW *r;
- - row_elt *e;
- -
- - if ( col < 0 || col >= A->n )
- - error(E_BOUNDS,"chase_col");
- - tmp_row = *row_num;
- - if ( tmp_row < 0 )
- - {
- - if ( A->start_row[col] > max_row )
- - {
- - tmp_row = -1;
- - tmp_idx = col;
- - return (row_elt *)NULL;
- - }
- - else
- - {
- - tmp_row = A->start_row[col];
- - tmp_idx = A->start_idx[col];
- - }
- - }
- - else
- - tmp_idx = *idx;
- -
- - old_row = tmp_row;
- - old_idx = tmp_idx;
- - while ( tmp_row >= 0 && tmp_row < max_row )
- - {
- - r = &(A->row[tmp_row]);
- - /* tmp_idx = sprow_idx2(r,col,tmp_idx); */
- - if ( tmp_idx < 0 || tmp_idx >= r->len ||
- - r->elt[tmp_idx].col != col )
- - {
- -#ifdef DEBUG
- - printf("chase_col:error: col = %d, row # = %d, idx = %d\n",
- - col, tmp_row, tmp_idx);
- - printf("chase_col:error: old_row = %d, old_idx = %d\n",
- - old_row, old_idx);
- - printf("chase_col:error: A =\n");
- - sp_dump(stdout,A);
- -#endif
- - error(E_INTERN,"chase_col");
- - }
- - e = &(r->elt[tmp_idx]);
- - old_row = tmp_row;
- - old_idx = tmp_idx;
- - tmp_row = e->nxt_row;
- - tmp_idx = e->nxt_idx;
- - }
- - if ( old_row > max_row )
- - {
- - old_row = -1;
- - old_idx = col;
- - e = (row_elt *)NULL;
- - }
- - else if ( tmp_row <= max_row && tmp_row >= 0 )
- - {
- - old_row = tmp_row;
- - old_idx = tmp_idx;
- - }
- -
- - *row_num = old_row;
- - if ( old_row >= 0 )
- - *idx = old_idx;
- - else
- - *idx = col;
- -
- - return e;
- -}
- -
- -/* chase_past -- as for chase_col except that we want the first
- - row whose row # >= min_row; -1 indicates no such row */
- -row_elt *chase_past(A, col, row_num, idx, min_row)
- -SPMAT *A;
- -int col, *row_num, *idx, min_row;
- -{
- - SPROW *r;
- - row_elt *e;
- - int tmp_idx, tmp_row;
- -
- - tmp_row = *row_num;
- - tmp_idx = *idx;
- - chase_col(A,col,&tmp_row,&tmp_idx,min_row);
- - if ( tmp_row < 0 ) /* use A->start_row[..] etc. */
- - {
- - if ( A->start_row[col] < 0 )
- - tmp_row = -1;
- - else
- - {
- - tmp_row = A->start_row[col];
- - tmp_idx = A->start_idx[col];
- - }
- - }
- - else if ( tmp_row < min_row )
- - {
- - r = &(A->row[tmp_row]);
- - if ( tmp_idx < 0 || tmp_idx >= r->len ||
- - r->elt[tmp_idx].col != col )
- - error(E_INTERN,"chase_past");
- - tmp_row = r->elt[tmp_idx].nxt_row;
- - tmp_idx = r->elt[tmp_idx].nxt_idx;
- - }
- -
- - *row_num = tmp_row;
- - *idx = tmp_idx;
- - if ( tmp_row < 0 )
- - e = (row_elt *)NULL;
- - else
- - {
- - if ( tmp_idx < 0 || tmp_idx >= A->row[tmp_row].len ||
- - A->row[tmp_row].elt[tmp_idx].col != col )
- - error(E_INTERN,"bump_col");
- - e = &(A->row[tmp_row].elt[tmp_idx]);
- - }
- -
- - return e;
- -}
- -
- -/* bump_col -- move along to next nonzero entry in column col after row_num
- - -- update row_num and idx */
- -row_elt *bump_col(A, col, row_num, idx)
- -SPMAT *A;
- -int col, *row_num, *idx;
- -{
- - SPROW *r;
- - row_elt *e;
- - int tmp_row, tmp_idx;
- -
- - tmp_row = *row_num;
- - tmp_idx = *idx;
- - /* printf("bump_col: col = %d, row# = %d, idx = %d\n",
- - col, *row_num, *idx); */
- - if ( tmp_row < 0 )
- - {
- - tmp_row = A->start_row[col];
- - tmp_idx = A->start_idx[col];
- - }
- - else
- - {
- - r = &(A->row[tmp_row]);
- - if ( tmp_idx < 0 || tmp_idx >= r->len ||
- - r->elt[tmp_idx].col != col )
- - error(E_INTERN,"bump_col");
- - e = &(r->elt[tmp_idx]);
- - tmp_row = e->nxt_row;
- - tmp_idx = e->nxt_idx;
- - }
- - if ( tmp_row < 0 )
- - {
- - e = (row_elt *)NULL;
- - tmp_idx = col;
- - }
- - else
- - {
- - if ( tmp_idx < 0 || tmp_idx >= A->row[tmp_row].len ||
- - A->row[tmp_row].elt[tmp_idx].col != col )
- - error(E_INTERN,"bump_col");
- - e = &(A->row[tmp_row].elt[tmp_idx]);
- - }
- - *row_num = tmp_row;
- - *idx = tmp_idx;
- -
- - return e;
- -}
- -
- -
- //GO.SYSIN DD spswap.c
- echo iter0.c 1>&2
- sed >iter0.c <<'//GO.SYSIN DD iter0.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -
- -/* iter0.c 14/09/93 */
- -
- -/* ITERATIVE METHODS - service functions */
- -
- -/* functions for creating and releasing ITER structures;
- - for memory information;
- - for getting some values from an ITER variable;
- - for changing values in an ITER variable;
- - see also iter.c
- -*/
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include "iter.h"
- -
- -
- -static char rcsid[] = "$Id: iter0.c,v 1.1 1994/01/13 05:38:23 des Exp $";
- -
- -
- -/* standard functions */
- -
- -/* standard information */
- -void iter_std_info(ip,nres,res,Bres)
- -ITER *ip;
- -double nres;
- -VEC *res, *Bres;
- -{
- - if (nres >= 0.0)
- - printf(" %d. residual = %g\n",ip->steps,nres);
- - else
- - printf(" %d. residual = %g (WARNING !!! should be >= 0) \n",
- - ip->steps,nres);
- -}
- -
- -/* standard stopping criterion */
- -int iter_std_stop_crit(ip, nres, res, Bres)
- -ITER *ip;
- -double nres;
- -VEC *res, *Bres;
- -{
- - /* save initial value of the residual in init_res */
- - if (ip->steps == 0)
- - ip->init_res = fabs(nres);
- -
- - /* standard stopping criterium */
- - if (nres <= ip->init_res*ip->eps) return TRUE;
- - return FALSE;
- -}
- -
- -
- -/* iter_get - create a new structure pointing to ITER */
- -
- -ITER *iter_get(lenb, lenx)
- -int lenb, lenx;
- -{
- - ITER *ip;
- -
- - if ((ip = NEW(ITER)) == (ITER *) NULL)
- - error(E_MEM,"iter_get");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_ITER,0,sizeof(ITER));
- - mem_numvar(TYPE_ITER,1);
- - }
- -
- - /* default values */
- -
- - ip->shared_x = FALSE;
- - ip->shared_b = FALSE;
- - ip->k = 0;
- - ip->limit = ITER_LIMIT_DEF;
- - ip->eps = ITER_EPS_DEF;
- - ip->steps = 0;
- -
- - if (lenb > 0) ip->b = v_get(lenb);
- - else ip->b = (VEC *)NULL;
- -
- - if (lenx > 0) ip->x = v_get(lenx);
- - else ip->x = (VEC *)NULL;
- -
- - ip->Ax = ip->A_par = NULL;
- - ip->ATx = ip->AT_par = NULL;
- - ip->Bx = ip->B_par = NULL;
- - ip->info = iter_std_info;
- - ip->stop_crit = iter_std_stop_crit;
- - ip->init_res = 0.0;
- -
- - return ip;
- -}
- -
- -
- -/* iter_free - release memory */
- -int iter_free(ip)
- -ITER *ip;
- -{
- - if (ip == (ITER *)NULL) return -1;
- -
- - if (mem_info_is_on()) {
- - mem_bytes(TYPE_ITER,sizeof(ITER),0);
- - mem_numvar(TYPE_ITER,-1);
- - }
- -
- - if ( !ip->shared_x && ip->x != NULL ) v_free(ip->x);
- - if ( !ip->shared_b && ip->b != NULL ) v_free(ip->b);
- -
- - free((char *)ip);
- -
- - return 0;
- -}
- -
- -ITER *iter_resize(ip,new_lenb,new_lenx)
- -ITER *ip;
- -int new_lenb, new_lenx;
- -{
- - VEC *old;
- -
- - if ( ip == (ITER *) NULL)
- - error(E_NULL,"iter_resize");
- -
- - old = ip->x;
- - ip->x = v_resize(ip->x,new_lenx);
- - if ( ip->shared_x && old != ip->x )
- - warning(WARN_SHARED_VEC,"iter_resize");
- - old = ip->b;
- - ip->b = v_resize(ip->b,new_lenb);
- - if ( ip->shared_b && old != ip->b )
- - warning(WARN_SHARED_VEC,"iter_resize");
- -
- - return ip;
- -}
- -
- -
- -/* print out ip structure - for diagnostic purposes mainly */
- -void iter_dump(fp,ip)
- -ITER *ip;
- -FILE *fp;
- -{
- - if (ip == NULL) {
- - fprintf(fp," ITER structure: NULL\n");
- - return;
- - }
- -
- - fprintf(fp,"\n ITER structure:\n");
- - fprintf(fp," ip->shared_x = %s, ip->shared_b = %s\n",
- - (ip->shared_x ? "TRUE" : "FALSE"),
- - (ip->shared_b ? "TRUE" : "FALSE") );
- - fprintf(fp," ip->k = %d, ip->limit = %d, ip->steps = %d, ip->eps = %g\n",
- - ip->k,ip->limit,ip->steps,ip->eps);
- - fprintf(fp," ip->x = 0x%p, ip->b = 0x%p\n",ip->x,ip->b);
- - fprintf(fp," ip->Ax = 0x%p, ip->A_par = 0x%p\n",ip->Ax,ip->A_par);
- - fprintf(fp," ip->ATx = 0x%p, ip->AT_par = 0x%p\n",ip->ATx,ip->AT_par);
- - fprintf(fp," ip->Bx = 0x%p, ip->B_par = 0x%p\n",ip->Bx,ip->B_par);
- - fprintf(fp," ip->info = 0x%p, ip->stop_crit = 0x%p, ip->init_res = %g\n",
- - ip->info,ip->stop_crit,ip->init_res);
- - fprintf(fp,"\n");
- -
- -}
- -
- -
- -/* copy the structure ip1 to ip2 preserving vectors x and b of ip2
- - (vectors x and b in ip2 are the same before and after iter_copy2)
- - if ip2 == NULL then a new structure is created with x and b being NULL
- - and other members are taken from ip1
- -*/
- -ITER *iter_copy2(ip1,ip2)
- -ITER *ip1, *ip2;
- -{
- - VEC *x, *b;
- - int shx, shb;
- -
- - if (ip1 == (ITER *)NULL)
- - error(E_NULL,"iter_copy2");
- -
- - if (ip2 == (ITER *)NULL) {
- - if ((ip2 = NEW(ITER)) == (ITER *) NULL)
- - error(E_MEM,"iter_copy2");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_ITER,0,sizeof(ITER));
- - mem_numvar(TYPE_ITER,1);
- - }
- - ip2->x = ip2->b = NULL;
- - ip2->shared_x = ip2->shared_x = FALSE;
- - }
- -
- - x = ip2->x;
- - b = ip2->b;
- - shb = ip2->shared_b;
- - shx = ip2->shared_x;
- - MEM_COPY(ip1,ip2,sizeof(ITER));
- - ip2->x = x;
- - ip2->b = b;
- - ip2->shared_x = shx;
- - ip2->shared_b = shb;
- -
- - return ip2;
- -}
- -
- -
- -/* copy the structure ip1 to ip2 copying also the vectors x and b */
- -ITER *iter_copy(ip1,ip2)
- -ITER *ip1, *ip2;
- -{
- - VEC *x, *b;
- -
- - if (ip1 == (ITER *)NULL)
- - error(E_NULL,"iter_copy");
- -
- - if (ip2 == (ITER *)NULL) {
- - if ((ip2 = NEW(ITER)) == (ITER *) NULL)
- - error(E_MEM,"iter_copy2");
- - else if (mem_info_is_on()) {
- - mem_bytes(TYPE_ITER,0,sizeof(ITER));
- - mem_numvar(TYPE_ITER,1);
- - }
- - }
- -
- - x = ip2->x;
- - b = ip2->b;
- -
- - MEM_COPY(ip1,ip2,sizeof(ITER));
- - if (ip1->x)
- - ip2->x = v_copy(ip1->x,x);
- - if (ip1->b)
- - ip2->b = v_copy(ip1->b,b);
- -
- - ip2->shared_x = ip2->shared_b = FALSE;
- -
- - return ip2;
- -}
- -
- -
- -/*** functions to generate sparse matrices with random entries ***/
- -
- -
- -/* iter_gen_sym -- generate symmetric positive definite
- - n x n matrix,
- - nrow - number of nonzero entries in a row
- - */
- -SPMAT *iter_gen_sym(n,nrow)
- -int n, nrow;
- -{
- - SPMAT *A;
- - VEC *u;
- - Real s1;
- - int i, j, k, k_max;
- -
- - if (nrow <= 1) nrow = 2;
- - /* nrow should be even */
- - if ((nrow & 1)) nrow -= 1;
- - A = sp_get(n,n,nrow);
- - u = v_get(A->m);
- - v_zero(u);
- - for ( i = 0; i < A->m; i++ )
- - {
- - k_max = ((rand() >> 8) % (nrow/2));
- - for ( k = 0; k <= k_max; k++ )
- - {
- - j = (rand() >> 8) % A->n;
- - s1 = mrand();
- - sp_set_val(A,i,j,s1);
- - sp_set_val(A,j,i,s1);
- - u->ve[i] += fabs(s1);
- - u->ve[j] += fabs(s1);
- - }
- - }
- - /* ensure that A is positive definite */
- - for ( i = 0; i < A->m; i++ )
- - sp_set_val(A,i,i,u->ve[i] + 1.0);
- -
- - V_FREE(u);
- - return A;
- -}
- -
- -
- -/* iter_gen_nonsym -- generate non-symmetric m x n sparse matrix, m >= n
- - nrow - number of entries in a row;
- - diag - number which is put in diagonal entries and then permuted
- - (if diag is zero then 1.0 is there)
- -*/
- -SPMAT *iter_gen_nonsym(m,n,nrow,diag)
- -int m, n, nrow;
- -double diag;
- -{
- - SPMAT *A;
- - PERM *px;
- - int i, j, k, k_max;
- - Real s1;
- -
- - if (nrow <= 1) nrow = 2;
- - if (diag == 0.0) diag = 1.0;
- - A = sp_get(m,n,nrow);
- - px = px_get(n);
- - for ( i = 0; i < A->m; i++ )
- - {
- - k_max = (rand() >> 8) % (nrow-1);
- - for ( k = 0; k <= k_max; k++ )
- - {
- - j = (rand() >> 8) % A->n;
- - s1 = mrand();
- - sp_set_val(A,i,j,-s1);
- - }
- - }
- - /* to make it likely that A is nonsingular, use pivot... */
- - for ( i = 0; i < 2*A->n; i++ )
- - {
- - j = (rand() >> 8) % A->n;
- - k = (rand() >> 8) % A->n;
- - px_transp(px,j,k);
- - }
- - for ( i = 0; i < A->n; i++ )
- - sp_set_val(A,i,px->pe[i],diag);
- -
- - PX_FREE(px);
- - return A;
- -}
- -
- -
- -/* iter_gen_nonsym -- generate non-symmetric positive definite
- - n x n sparse matrix;
- - nrow - number of entries in a row
- -*/
- -SPMAT *iter_gen_nonsym_posdef(n,nrow)
- -int n, nrow;
- -{
- - SPMAT *A;
- - PERM *px;
- - VEC *u;
- - int i, j, k, k_max;
- - Real s1;
- -
- - if (nrow <= 1) nrow = 2;
- - A = sp_get(n,n,nrow);
- - px = px_get(n);
- - u = v_get(A->m);
- - v_zero(u);
- - for ( i = 0; i < A->m; i++ )
- - {
- - k_max = (rand() >> 8) % (nrow-1);
- - for ( k = 0; k <= k_max; k++ )
- - {
- - j = (rand() >> 8) % A->n;
- - s1 = mrand();
- - sp_set_val(A,i,j,-s1);
- - u->ve[i] += fabs(s1);
- - }
- - }
- - /* ensure that A is positive definite */
- - for ( i = 0; i < A->m; i++ )
- - sp_set_val(A,i,i,u->ve[i] + 1.0);
- -
- - PX_FREE(px);
- - V_FREE(u);
- - return A;
- -}
- -
- -
- //GO.SYSIN DD iter0.c
- echo itersym.c 1>&2
- sed >itersym.c <<'//GO.SYSIN DD itersym.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -
- -/* itersym.c 17/09/93 */
- -
- -
- -/*
- - ITERATIVE METHODS - implementation of several iterative methods;
- - see also iter0.c
- - */
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include "matrix.h"
- -#include "matrix2.h"
- -#include "sparse.h"
- -#include "iter.h"
- -
- -static char rcsid[] = "$Id: itersym.c,v 1.1 1994/01/13 05:38:59 des Exp $";
- -
- -
- -#ifdef ANSI_C
- -VEC *spCHsolve(SPMAT *,VEC *,VEC *);
- -VEC *trieig(VEC *,VEC *,MAT *);
- -#else
- -VEC *spCHsolve();
- -VEC *trieig();
- -#endif
- -
- -
- -
- -/* iter_spcg -- a simple interface to iter_cg() which uses sparse matrix
- - data structures
- - -- assumes that LLT contains the Cholesky factorisation of the
- - actual preconditioner;
- - use always as follows:
- - x = iter_spcg(A,LLT,b,eps,x,limit,steps);
- - or
- - x = iter_spcg(A,LLT,b,eps,VNULL,limit,steps);
- - In the second case the solution vector is created.
- - */
- -VEC *iter_spcg(A,LLT,b,eps,x,limit,steps)
- -SPMAT *A, *LLT;
- -VEC *b, *x;
- -double eps;
- -int *steps, limit;
- -{
- - ITER *ip;
- -
- - ip = iter_get(0,0);
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *)A;
- - ip->Bx = (Fun_Ax) spCHsolve;
- - ip->B_par = (void *)LLT;
- - ip->info = (Fun_info) NULL;
- - ip->b = b;
- - ip->eps = eps;
- - ip->limit = limit;
- - ip->x = x;
- - iter_cg(ip);
- - x = ip->x;
- - if (steps) *steps = ip->steps;
- - ip->shared_x = ip->shared_b = TRUE;
- - iter_free(ip); /* release only ITER structure */
- - return x;
- -}
- -
- -/*
- - Conjugate gradients method;
- - */
- -VEC *iter_cg(ip)
- -ITER *ip;
- -{
- - static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
- - Real alpha, beta, inner, old_inner, nres;
- - VEC *rr; /* rr == r or rr == z */
- -
- - if (ip == INULL)
- - error(E_NULL,"iter_cg");
- - if (!ip->Ax || !ip->b)
- - error(E_NULL,"iter_cg");
- - if ( ip->x == ip->b )
- - error(E_INSITU,"iter_cg");
- - if (!ip->stop_crit)
- - error(E_NULL,"iter_cg");
- -
- - if ( ip->eps <= 0.0 )
- - ip->eps = MACHEPS;
- -
- - r = v_resize(r,ip->b->dim);
- - p = v_resize(p,ip->b->dim);
- - q = v_resize(q,ip->b->dim);
- -
- - MEM_STAT_REG(r,TYPE_VEC);
- - MEM_STAT_REG(p,TYPE_VEC);
- - MEM_STAT_REG(q,TYPE_VEC);
- -
- - if (ip->Bx != (Fun_Ax)NULL) {
- - z = v_resize(z,ip->b->dim);
- - MEM_STAT_REG(z,TYPE_VEC);
- - rr = z;
- - }
- - else rr = r;
- -
- - if (ip->x != VNULL) {
- - if (ip->x->dim != ip->b->dim)
- - error(E_SIZES,"iter_cg");
- - ip->Ax(ip->A_par,ip->x,p); /* p = A*x */
- - v_sub(ip->b,p,r); /* r = b - A*x */
- - }
- - else { /* ip->x == 0 */
- - ip->x = v_get(ip->b->dim);
- - ip->shared_x = FALSE;
- - v_copy(ip->b,r);
- - }
- -
- - old_inner = 0.0;
- - for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
- - {
- - if ( ip->Bx )
- - (ip->Bx)(ip->B_par,r,rr); /* rr = B*r */
- -
- - inner = in_prod(rr,r);
- - nres = sqrt(fabs(inner));
- - if (ip->info) ip->info(ip,nres,r,rr);
- - if ( ip->stop_crit(ip,nres,r,rr) ) break;
- -
- - if ( ip->steps ) /* if ( ip->steps > 0 ) ... */
- - {
- - beta = inner/old_inner;
- - p = v_mltadd(rr,p,beta,p);
- - }
- - else /* if ( ip->steps == 0 ) ... */
- - {
- - beta = 0.0;
- - p = v_copy(rr,p);
- - old_inner = 0.0;
- - }
- - (ip->Ax)(ip->A_par,p,q); /* q = A*p */
- - alpha = inner/in_prod(p,q);
- - v_mltadd(ip->x,p,alpha,ip->x);
- - v_mltadd(r,q,-alpha,r);
- - old_inner = inner;
- - }
- -
- - return ip->x;
- -}
- -
- -
- -
- -/* iter_lanczos -- raw lanczos algorithm -- no re-orthogonalisation
- - -- creates T matrix of size == m,
- - but no larger than before beta_k == 0
- - -- uses passed routine to do matrix-vector multiplies */
- -void iter_lanczos(ip,a,b,beta2,Q)
- -ITER *ip;
- -VEC *a, *b;
- -Real *beta2;
- -MAT *Q;
- -{
- - int j;
- - static VEC *v = VNULL, *w = VNULL, *tmp = VNULL;
- - Real alpha, beta, c;
- -
- - if ( ! ip )
- - error(E_NULL,"iter_lanczos");
- - if ( ! ip->Ax || ! ip->x || ! a || ! b )
- - error(E_NULL,"iter_lanczos");
- - if ( ip->k <= 0 )
- - error(E_BOUNDS,"iter_lanczos");
- - if ( Q && ( Q->n < ip->x->dim || Q->m < ip->k ) )
- - error(E_SIZES,"iter_lanczos");
- -
- - a = v_resize(a,(u_int)ip->k);
- - b = v_resize(b,(u_int)(ip->k-1));
- - v = v_resize(v,ip->x->dim);
- - w = v_resize(w,ip->x->dim);
- - tmp = v_resize(tmp,ip->x->dim);
- - MEM_STAT_REG(v,TYPE_VEC);
- - MEM_STAT_REG(w,TYPE_VEC);
- - MEM_STAT_REG(tmp,TYPE_VEC);
- -
- - beta = 1.0;
- - v_zero(a);
- - v_zero(b);
- - if (Q) m_zero(Q);
- -
- - /* normalise x as w */
- - c = v_norm2(ip->x);
- - if (c <= MACHEPS) { /* ip->x == 0 */
- - *beta2 = 0.0;
- - return;
- - }
- - else
- - sv_mlt(1.0/c,ip->x,w);
- -
- - (ip->Ax)(ip->A_par,w,v);
- -
- - for ( j = 0; j < ip->k; j++ )
- - {
- - /* store w in Q if Q not NULL */
- - if ( Q ) set_row(Q,j,w);
- -
- - alpha = in_prod(w,v);
- - a->ve[j] = alpha;
- - v_mltadd(v,w,-alpha,v);
- - beta = v_norm2(v);
- - if ( beta == 0.0 )
- - {
- - *beta2 = 0.0;
- - return;
- - }
- -
- - if ( j < ip->k-1 )
- - b->ve[j] = beta;
- - v_copy(w,tmp);
- - sv_mlt(1/beta,v,w);
- - sv_mlt(-beta,tmp,v);
- - (ip->Ax)(ip->A_par,w,tmp);
- - v_add(v,tmp,v);
- - }
- - *beta2 = beta;
- -
- -}
- -
- -/* iter_splanczos -- version that uses sparse matrix data structure */
- -void iter_splanczos(A,m,x0,a,b,beta2,Q)
- -SPMAT *A;
- -int m;
- -VEC *x0, *a, *b;
- -Real *beta2;
- -MAT *Q;
- -{
- - ITER *ip;
- -
- - ip = iter_get(0,0);
- - ip->shared_x = ip->shared_b = TRUE;
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *) A;
- - ip->x = x0;
- - ip->k = m;
- - iter_lanczos(ip,a,b,beta2,Q);
- - iter_free(ip); /* release only ITER structure */
- -}
- -
- -
- -
- -extern double frexp(), ldexp();
- -
- -/* product -- returns the product of a long list of numbers
- - -- answer stored in mant (mantissa) and expt (exponent) */
- -static double product(a,offset,expt)
- -VEC *a;
- -double offset;
- -int *expt;
- -{
- - Real mant, tmp_fctr;
- - int i, tmp_expt;
- -
- - if ( ! a )
- - error(E_NULL,"product");
- -
- - mant = 1.0;
- - *expt = 0;
- - if ( offset == 0.0 )
- - for ( i = 0; i < a->dim; i++ )
- - {
- - mant *= frexp(a->ve[i],&tmp_expt);
- - *expt += tmp_expt;
- - if ( ! (i % 10) )
- - {
- - mant = frexp(mant,&tmp_expt);
- - *expt += tmp_expt;
- - }
- - }
- - else
- - for ( i = 0; i < a->dim; i++ )
- - {
- - tmp_fctr = a->ve[i] - offset;
- - tmp_fctr += (tmp_fctr > 0.0 ) ? -MACHEPS*offset :
- - MACHEPS*offset;
- - mant *= frexp(tmp_fctr,&tmp_expt);
- - *expt += tmp_expt;
- - if ( ! (i % 10) )
- - {
- - mant = frexp(mant,&tmp_expt);
- - *expt += tmp_expt;
- - }
- - }
- -
- - mant = frexp(mant,&tmp_expt);
- - *expt += tmp_expt;
- -
- - return mant;
- -}
- -
- -/* product2 -- returns the product of a long list of numbers
- - -- answer stored in mant (mantissa) and expt (exponent) */
- -static double product2(a,k,expt)
- -VEC *a;
- -int k; /* entry of a to leave out */
- -int *expt;
- -{
- - Real mant, mu, tmp_fctr;
- - int i, tmp_expt;
- -
- - if ( ! a )
- - error(E_NULL,"product2");
- - if ( k < 0 || k >= a->dim )
- - error(E_BOUNDS,"product2");
- -
- - mant = 1.0;
- - *expt = 0;
- - mu = a->ve[k];
- - for ( i = 0; i < a->dim; i++ )
- - {
- - if ( i == k )
- - continue;
- - tmp_fctr = a->ve[i] - mu;
- - tmp_fctr += ( tmp_fctr > 0.0 ) ? -MACHEPS*mu : MACHEPS*mu;
- - mant *= frexp(tmp_fctr,&tmp_expt);
- - *expt += tmp_expt;
- - if ( ! (i % 10) )
- - {
- - mant = frexp(mant,&tmp_expt);
- - *expt += tmp_expt;
- - }
- - }
- - mant = frexp(mant,&tmp_expt);
- - *expt += tmp_expt;
- -
- - return mant;
- -}
- -
- -/* dbl_cmp -- comparison function to pass to qsort() */
- -static int dbl_cmp(x,y)
- -Real *x, *y;
- -{
- - Real tmp;
- -
- - tmp = *x - *y;
- - return (tmp > 0 ? 1 : tmp < 0 ? -1: 0);
- -}
- -
- -/* iter_lanczos2 -- lanczos + error estimate for every e-val
- - -- uses Cullum & Willoughby approach, Sparse Matrix Proc. 1978
- - -- returns multiple e-vals where multiple e-vals may not exist
- - -- returns evals vector */
- -VEC *iter_lanczos2(ip,evals,err_est)
- -ITER *ip; /* ITER structure */
- -VEC *evals; /* eigenvalue vector */
- -VEC *err_est; /* error estimates of eigenvalues */
- -{
- - VEC *a;
- - static VEC *b=VNULL, *a2=VNULL, *b2=VNULL;
- - Real beta, pb_mant, det_mant, det_mant1, det_mant2;
- - int i, pb_expt, det_expt, det_expt1, det_expt2;
- -
- - if ( ! ip )
- - error(E_NULL,"iter_lanczos2");
- - if ( ! ip->Ax || ! ip->x )
- - error(E_NULL,"iter_lanczos2");
- - if ( ip->k <= 0 )
- - error(E_RANGE,"iter_lanczos2");
- -
- - a = evals;
- - a = v_resize(a,(u_int)ip->k);
- - b = v_resize(b,(u_int)(ip->k-1));
- - MEM_STAT_REG(b,TYPE_VEC);
- -
- - iter_lanczos(ip,a,b,&beta,MNULL);
- -
- - /* printf("# beta =%g\n",beta); */
- - pb_mant = 0.0;
- - if ( err_est )
- - {
- - pb_mant = product(b,(double)0.0,&pb_expt);
- - /* printf("# pb_mant = %g, pb_expt = %d\n",pb_mant, pb_expt); */
- - }
- -
- - /* printf("# diags =\n"); v_output(a); */
- - /* printf("# off diags =\n"); v_output(b); */
- - a2 = v_resize(a2,a->dim - 1);
- - b2 = v_resize(b2,b->dim - 1);
- - MEM_STAT_REG(a2,TYPE_VEC);
- - MEM_STAT_REG(b2,TYPE_VEC);
- - for ( i = 0; i < a2->dim - 1; i++ )
- - {
- - a2->ve[i] = a->ve[i+1];
- - b2->ve[i] = b->ve[i+1];
- - }
- - a2->ve[a2->dim-1] = a->ve[a2->dim];
- -
- - trieig(a,b,MNULL);
- -
- - /* sort evals as a courtesy */
- - qsort((void *)(a->ve),(int)(a->dim),sizeof(Real),(int (*)())dbl_cmp);
- -
- - /* error estimates */
- - if ( err_est )
- - {
- - err_est = v_resize(err_est,(u_int)ip->k);
- -
- - trieig(a2,b2,MNULL);
- - /* printf("# a =\n"); v_output(a); */
- - /* printf("# a2 =\n"); v_output(a2); */
- -
- - for ( i = 0; i < a->dim; i++ )
- - {
- - det_mant1 = product2(a,i,&det_expt1);
- - det_mant2 = product(a2,(double)a->ve[i],&det_expt2);
- - /* printf("# det_mant1=%g, det_expt1=%d\n",
- - det_mant1,det_expt1); */
- - /* printf("# det_mant2=%g, det_expt2=%d\n",
- - det_mant2,det_expt2); */
- - if ( det_mant1 == 0.0 )
- - { /* multiple e-val of T */
- - err_est->ve[i] = 0.0;
- - continue;
- - }
- - else if ( det_mant2 == 0.0 )
- - {
- - err_est->ve[i] = HUGE;
- - continue;
- - }
- - if ( (det_expt1 + det_expt2) % 2 )
- - /* if odd... */
- - det_mant = sqrt(2.0*fabs(det_mant1*det_mant2));
- - else /* if even... */
- - det_mant = sqrt(fabs(det_mant1*det_mant2));
- - det_expt = (det_expt1+det_expt2)/2;
- - err_est->ve[i] = fabs(beta*
- - ldexp(pb_mant/det_mant,pb_expt-det_expt));
- - }
- - }
- -
- - return a;
- -}
- -
- -/* iter_splanczos2 -- version of iter_lanczos2() that uses sparse matrix data
- - structure */
- -
- -VEC *iter_splanczos2(A,m,x0,evals,err_est)
- -SPMAT *A;
- -int m;
- -VEC *x0; /* initial vector */
- -VEC *evals; /* eigenvalue vector */
- -VEC *err_est; /* error estimates of eigenvalues */
- -{
- - ITER *ip;
- - VEC *a;
- -
- - ip = iter_get(0,0);
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *) A;
- - ip->x = x0;
- - ip->k = m;
- - a = iter_lanczos2(ip,evals,err_est);
- - ip->shared_x = ip->shared_b = TRUE;
- - iter_free(ip); /* release only ITER structure */
- - return a;
- -}
- -
- -
- -
- -
- -/*
- - Conjugate gradient method
- - Another variant - mainly for testing
- - */
- -
- -VEC *iter_cg1(ip)
- -ITER *ip;
- -{
- - static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
- - Real alpha;
- - double inner,nres;
- - VEC *rr; /* rr == r or rr == z */
- -
- - if (ip == INULL)
- - error(E_NULL,"iter_cg");
- - if (!ip->Ax || !ip->b)
- - error(E_NULL,"iter_cg");
- - if ( ip->x == ip->b )
- - error(E_INSITU,"iter_cg");
- - if (!ip->stop_crit)
- - error(E_NULL,"iter_cg");
- -
- - if ( ip->eps <= 0.0 )
- - ip->eps = MACHEPS;
- -
- - r = v_resize(r,ip->b->dim);
- - p = v_resize(p,ip->b->dim);
- - q = v_resize(q,ip->b->dim);
- -
- - MEM_STAT_REG(r,TYPE_VEC);
- - MEM_STAT_REG(p,TYPE_VEC);
- - MEM_STAT_REG(q,TYPE_VEC);
- -
- - if (ip->Bx != (Fun_Ax)NULL) {
- - z = v_resize(z,ip->b->dim);
- - MEM_STAT_REG(z,TYPE_VEC);
- - rr = z;
- - }
- - else rr = r;
- -
- - if (ip->x != VNULL) {
- - if (ip->x->dim != ip->b->dim)
- - error(E_SIZES,"iter_cg");
- - ip->Ax(ip->A_par,ip->x,p); /* p = A*x */
- - v_sub(ip->b,p,r); /* r = b - A*x */
- - }
- - else { /* ip->x == 0 */
- - ip->x = v_get(ip->b->dim);
- - ip->shared_x = FALSE;
- - v_copy(ip->b,r);
- - }
- -
- - if (ip->Bx) (ip->Bx)(ip->B_par,r,p);
- - else v_copy(r,p);
- -
- - inner = in_prod(p,r);
- - nres = sqrt(fabs(inner));
- - if (ip->info) ip->info(ip,nres,r,p);
- - if ( ip->stop_crit(ip,nres,r,p) ) return ip->x;
- -
- - for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
- - {
- - ip->Ax(ip->A_par,p,q);
- - inner = in_prod(q,p);
- - if (inner <= 0.0) {
- - warning(WARN_RES_LESS_0,"iter_cg");
- - break;
- - }
- - alpha = in_prod(p,r)/inner;
- - v_mltadd(ip->x,p,alpha,ip->x);
- - v_mltadd(r,q,-alpha,r);
- -
- - rr = r;
- - if (ip->Bx) {
- - ip->Bx(ip->B_par,r,z);
- - rr = z;
- - }
- -
- - nres = in_prod(r,rr);
- - if (nres < 0.0) {
- - warning(WARN_RES_LESS_0,"iter_cg");
- - break;
- - }
- - nres = sqrt(fabs(nres));
- - if (ip->info) ip->info(ip,nres,r,z);
- - if ( ip->stop_crit(ip,nres,r,z) ) break;
- -
- - alpha = -in_prod(rr,q)/inner;
- - v_mltadd(rr,p,alpha,p);
- -
- - }
- -
- - return ip->x;
- -}
- -
- //GO.SYSIN DD itersym.c
- echo iternsym.c 1>&2
- sed >iternsym.c <<'//GO.SYSIN DD iternsym.c' 's/^-//'
- -
- -/**************************************************************************
- -**
- -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
- -**
- -** Meschach Library
- -**
- -** This Meschach Library is provided "as is" without any express
- -** or implied warranty of any kind with respect to this software.
- -** In particular the authors shall not be liable for any direct,
- -** indirect, special, incidental or consequential damages arising
- -** in any way from use of the software.
- -**
- -** Everyone is granted permission to copy, modify and redistribute this
- -** Meschach Library, provided:
- -** 1. All copies contain this copyright notice.
- -** 2. All modified copies shall carry a notice stating who
- -** made the last modification and the date of such modification.
- -** 3. No charge is made for this software or works derived from it.
- -** This clause shall not be construed as constraining other software
- -** distributed on the same medium as this software, nor is a
- -** distribution fee considered a charge.
- -**
- -***************************************************************************/
- -
- -
- -/* iter.c 17/09/93 */
- -
- -/*
- - ITERATIVE METHODS - implementation of several iterative methods;
- - see also iter0.c
- -*/
- -
- -#include <stdio.h>
- -#include <math.h>
- -#include "matrix.h"
- -#include "matrix2.h"
- -#include "sparse.h"
- -#include "iter.h"
- -
- -static char rcsid[] = "$Id: iternsym.c,v 1.2 1994/02/02 06:02:20 des Exp $";
- -
- -
- -#ifdef ANSI_C
- -VEC *spCHsolve(SPMAT *,VEC *,VEC *);
- -#else
- -VEC *spCHsolve();
- -#endif
- -
- -
- -/*
- - iter_cgs -- uses CGS to compute a solution x to A.x=b
- -*/
- -
- -VEC *iter_cgs(ip,r0)
- -ITER *ip;
- -VEC *r0;
- -{
- - static VEC *p = VNULL, *q = VNULL, *r = VNULL, *u = VNULL;
- - static VEC *v = VNULL, *z = VNULL;
- - VEC *tmp;
- - Real alpha, beta, nres, rho, old_rho, sigma, inner;
- -
- - if (ip == INULL)
- - error(E_NULL,"iter_cgs");
- - if (!ip->Ax || !ip->b || !r0)
- - error(E_NULL,"iter_cgs");
- - if ( ip->x == ip->b )
- - error(E_INSITU,"iter_cgs");
- - if (!ip->stop_crit)
- - error(E_NULL,"iter_cgs");
- - if ( r0->dim != ip->b->dim )
- - error(E_SIZES,"iter_cgs");
- -
- - if ( ip->eps <= 0.0 )
- - ip->eps = MACHEPS;
- -
- - p = v_resize(p,ip->b->dim);
- - q = v_resize(q,ip->b->dim);
- - r = v_resize(r,ip->b->dim);
- - u = v_resize(u,ip->b->dim);
- - v = v_resize(v,ip->b->dim);
- -
- - MEM_STAT_REG(p,TYPE_VEC);
- - MEM_STAT_REG(q,TYPE_VEC);
- - MEM_STAT_REG(r,TYPE_VEC);
- - MEM_STAT_REG(u,TYPE_VEC);
- - MEM_STAT_REG(v,TYPE_VEC);
- -
- - if (ip->Bx) {
- - z = v_resize(z,ip->b->dim);
- - MEM_STAT_REG(z,TYPE_VEC);
- - }
- -
- - if (ip->x != VNULL) {
- - if (ip->x->dim != ip->b->dim)
- - error(E_SIZES,"iter_cgs");
- - ip->Ax(ip->A_par,ip->x,v); /* v = A*x */
- - if (ip->Bx) {
- - v_sub(ip->b,v,v); /* v = b - A*x */
- - (ip->Bx)(ip->B_par,v,r); /* r = B*(b-A*x) */
- - }
- - else v_sub(ip->b,v,r); /* r = b-A*x */
- - }
- - else { /* ip->x == 0 */
- - ip->x = v_get(ip->b->dim); /* x == 0 */
- - ip->shared_x = FALSE;
- - if (ip->Bx) (ip->Bx)(ip->B_par,ip->b,r); /* r = B*b */
- - else v_copy(ip->b,r); /* r = b */
- - }
- -
- - v_zero(p);
- - v_zero(q);
- - old_rho = 1.0;
- -
- - for (ip->steps = 0; ip->steps <= ip->limit; ip->steps++) {
- -
- - inner = in_prod(r,r);
- - nres = sqrt(fabs(inner));
- -
- - if (ip->info) ip->info(ip,nres,r,VNULL);
- - if ( ip->stop_crit(ip,nres,r,VNULL) ) break;
- -
- - rho = in_prod(r0,r);
- - if ( old_rho == 0.0 )
- - error(E_SING,"iter_cgs");
- - beta = rho/old_rho;
- - v_mltadd(r,q,beta,u);
- - v_mltadd(q,p,beta,v);
- - v_mltadd(u,v,beta,p);
- -
- - (ip->Ax)(ip->A_par,p,q);
- - if (ip->Bx) {
- - (ip->Bx)(ip->B_par,q,z);
- - tmp = z;
- - }
- - else tmp = q;
- -
- - sigma = in_prod(r0,tmp);
- - if ( sigma == 0.0 )
- - error(E_SING,"iter_cgs");
- - alpha = rho/sigma;
- - v_mltadd(u,tmp,-alpha,q);
- - v_add(u,q,v);
- -
- - (ip->Ax)(ip->A_par,v,u);
- - if (ip->Bx) {
- - (ip->Bx)(ip->B_par,u,z);
- - tmp = z;
- - }
- - else tmp = u;
- -
- - v_mltadd(r,tmp,-alpha,r);
- - v_mltadd(ip->x,v,alpha,ip->x);
- -
- - old_rho = rho;
- - }
- -
- - return ip->x;
- -}
- -
- -
- -
- -/* iter_spcgs -- simple interface for SPMAT data structures
- - use always as follows:
- - x = iter_spcgs(A,B,b,r0,tol,x,limit,steps);
- - or
- - x = iter_spcgs(A,B,b,r0,tol,VNULL,limit,steps);
- - In the second case the solution vector is created.
- - If B is not NULL then it is a preconditioner.
- -*/
- -VEC *iter_spcgs(A,B,b,r0,tol,x,limit,steps)
- -SPMAT *A, *B;
- -VEC *b, *r0, *x;
- -double tol;
- -int *steps,limit;
- -{
- - ITER *ip;
- -
- - ip = iter_get(0,0);
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *) A;
- - if (B) {
- - ip->Bx = (Fun_Ax) sp_mv_mlt;
- - ip->B_par = (void *) B;
- - }
- - else {
- - ip->Bx = (Fun_Ax) NULL;
- - ip->B_par = NULL;
- - }
- - ip->info = (Fun_info) NULL;
- - ip->limit = limit;
- - ip->b = b;
- - ip->eps = tol;
- - ip->x = x;
- - iter_cgs(ip,r0);
- - x = ip->x;
- - if (steps) *steps = ip->steps;
- - ip->shared_x = ip->shared_b = TRUE;
- - iter_free(ip); /* release only ITER structure */
- - return x;
- -
- -}
- -
- -/*
- - Routine for performing LSQR -- the least squares QR algorithm
- - of Paige and Saunders:
- - "LSQR: an algorithm for sparse linear equations and
- - sparse least squares", ACM Trans. Math. Soft., v. 8
- - pp. 43--71 (1982)
- - */
- -/* lsqr -- sparse CG-like least squares routine:
- - -- finds min_x ||A.x-b||_2 using A defined through A & AT
- - -- returns x (if x != NULL) */
- -VEC *iter_lsqr(ip)
- -ITER *ip;
- -{
- - static VEC *u = VNULL, *v = VNULL, *w = VNULL, *tmp = VNULL;
- - Real alpha, beta, phi, phi_bar;
- - Real rho, rho_bar, rho_max, theta, nres;
- - Real s, c; /* for Givens' rotations */
- - int m, n;
- -
- - if ( ! ip || ! ip->b || !ip->Ax || !ip->ATx )
- - error(E_NULL,"iter_lsqr");
- - if ( ip->x == ip->b )
- - error(E_INSITU,"iter_lsqr");
- - if (!ip->stop_crit || !ip->x)
- - error(E_NULL,"iter_lsqr");
- -
- - if ( ip->eps <= 0.0 )
- - ip->eps = MACHEPS;
- -
- - m = ip->b->dim;
- - n = ip->x->dim;
- -
- - u = v_resize(u,(u_int)m);
- - v = v_resize(v,(u_int)n);
- - w = v_resize(w,(u_int)n);
- - tmp = v_resize(tmp,(u_int)n);
- -
- - MEM_STAT_REG(u,TYPE_VEC);
- - MEM_STAT_REG(v,TYPE_VEC);
- - MEM_STAT_REG(w,TYPE_VEC);
- - MEM_STAT_REG(tmp,TYPE_VEC);
- -
- - if (ip->x != VNULL) {
- - ip->Ax(ip->A_par,ip->x,u); /* u = A*x */
- - v_sub(ip->b,u,u); /* u = b-A*x */
- - }
- - else { /* ip->x == 0 */
- - ip->x = v_get(ip->b->dim);
- - ip->shared_x = FALSE;
- - v_copy(ip->b,u); /* u = b */
- - }
- -
- - beta = v_norm2(u);
- - if ( beta == 0.0 )
- - return ip->x;
- - sv_mlt(1.0/beta,u,u);
- - (ip->ATx)(ip->AT_par,u,v);
- - alpha = v_norm2(v);
- - if ( alpha == 0.0 )
- - return ip->x;
- - sv_mlt(1.0/alpha,v,v);
- - v_copy(v,w);
- - phi_bar = beta;
- - rho_bar = alpha;
- -
- - rho_max = 1.0;
- - for (ip->steps = 0; ip->steps <= ip->limit; ip->steps++) {
- -
- - tmp = v_resize(tmp,m);
- - (ip->Ax)(ip->A_par,v,tmp);
- -
- - v_mltadd(tmp,u,-alpha,u);
- - beta = v_norm2(u);
- - sv_mlt(1.0/beta,u,u);
- -
- - tmp = v_resize(tmp,n);
- - (ip->ATx)(ip->AT_par,u,tmp);
- - v_mltadd(tmp,v,-beta,v);
- - alpha = v_norm2(v);
- - sv_mlt(1.0/alpha,v,v);
- -
- - rho = sqrt(rho_bar*rho_bar+beta*beta);
- - if ( rho > rho_max )
- - rho_max = rho;
- - c = rho_bar/rho;
- - s = beta/rho;
- - theta = s*alpha;
- - rho_bar = -c*alpha;
- - phi = c*phi_bar;
- - phi_bar = s*phi_bar;
- -
- - /* update ip->x & w */
- - if ( rho == 0.0 )
- - error(E_SING,"iter_lsqr");
- - v_mltadd(ip->x,w,phi/rho,ip->x);
- - v_mltadd(v,w,-theta/rho,w);
- -
- - nres = fabs(phi_bar*alpha*c)*rho_max;
- -
- - if (ip->info) ip->info(ip,nres,w,VNULL);
- - if ( ip->stop_crit(ip,nres,w,VNULL) ) break;
- - }
- -
- - return ip->x;
- -}
- -
- -/* iter_splsqr -- simple interface for SPMAT data structures */
- -VEC *iter_splsqr(A,b,tol,x,limit,steps)
- -SPMAT *A;
- -VEC *b, *x;
- -double tol;
- -int *steps,limit;
- -{
- - ITER *ip;
- -
- - ip = iter_get(0,0);
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *) A;
- - ip->ATx = (Fun_Ax) sp_vm_mlt;
- - ip->AT_par = (void *) A;
- - ip->Bx = (Fun_Ax) NULL;
- - ip->B_par = NULL;
- -
- - ip->info = (Fun_info) NULL;
- - ip->limit = limit;
- - ip->b = b;
- - ip->eps = tol;
- - ip->x = x;
- - iter_lsqr(ip);
- - x = ip->x;
- - if (steps) *steps = ip->steps;
- - ip->shared_x = ip->shared_b = TRUE;
- - iter_free(ip); /* release only ITER structure */
- - return x;
- -}
- -
- -
- -
- -/* iter_arnoldi -- an implementation of the Arnoldi method;
- - iterative refinement is applied.
- -*/
- -MAT *iter_arnoldi_iref(ip,h_rem,Q,H)
- -ITER *ip;
- -Real *h_rem;
- -MAT *Q, *H;
- -{
- - static VEC *u=VNULL, *r=VNULL, *s=VNULL, *tmp=VNULL;
- - VEC v; /* auxiliary vector */
- - int i,j;
- - Real h_val, c;
- -
- - if (ip == INULL)
- - error(E_NULL,"iter_arnoldi_iref");
- - if ( ! ip->Ax || ! Q || ! ip->x )
- - error(E_NULL,"iter_arnoldi_iref");
- - if ( ip->k <= 0 )
- - error(E_BOUNDS,"iter_arnoldi_iref");
- - if ( Q->n != ip->x->dim || Q->m != ip->k )
- - error(E_SIZES,"iter_arnoldi_iref");
- -
- - m_zero(Q);
- - H = m_resize(H,ip->k,ip->k);
- - m_zero(H);
- -
- - u = v_resize(u,ip->x->dim);
- - tmp = v_resize(tmp,ip->x->dim);
- - r = v_resize(r,ip->k);
- - s = v_resize(s,ip->k);
- - MEM_STAT_REG(u,TYPE_VEC);
- - MEM_STAT_REG(tmp,TYPE_VEC);
- - MEM_STAT_REG(r,TYPE_VEC);
- - MEM_STAT_REG(s,TYPE_VEC);
- -
- - v.dim = v.max_dim = ip->x->dim;
- -
- - c = v_norm2(ip->x);
- - if ( c <= 0.0)
- - return H;
- - else {
- - v.ve = Q->me[0];
- - sv_mlt(1.0/c,ip->x,&v);
- - }
- -
- - v_zero(r);
- - v_zero(s);
- - for ( i = 0; i < ip->k; i++ )
- - {
- - v.ve = Q->me[i];
- - u = (ip->Ax)(ip->A_par,&v,u);
- - v_zero(tmp);
- - for (j = 0; j <= i; j++) {
- - v.ve = Q->me[j];
- - r->ve[j] = in_prod(&v,u);
- - v_mltadd(tmp,&v,r->ve[j],tmp);
- - }
- - v_sub(u,tmp,u);
- - h_val = v_norm2(u);
- - /* if u == 0 then we have an exact subspace */
- - if ( h_val <= 0.0 )
- - {
- - *h_rem = h_val;
- - return H;
- - }
- - /* iterative refinement -- ensures near orthogonality */
- - do {
- - v_zero(tmp);
- - for (j = 0; j <= i; j++) {
- - v.ve = Q->me[j];
- - s->ve[j] = in_prod(&v,u);
- - v_mltadd(tmp,&v,s->ve[j],tmp);
- - }
- - v_sub(u,tmp,u);
- - v_add(r,s,r);
- - } while ( v_norm2(s) > 0.1*(h_val = v_norm2(u)) );
- - /* now that u is nearly orthogonal to Q, update H */
- - set_col(H,i,r);
- - /* check once again if h_val is zero */
- - if ( h_val <= 0.0 )
- - {
- - *h_rem = h_val;
- - return H;
- - }
- - if ( i == ip->k-1 )
- - {
- - *h_rem = h_val;
- - continue;
- - }
- - /* H->me[i+1][i] = h_val; */
- - m_set_val(H,i+1,i,h_val);
- - v.ve = Q->me[i+1];
- - sv_mlt(1.0/h_val,u,&v);
- - }
- -
- - return H;
- -}
- -
- -/* iter_arnoldi -- an implementation of the Arnoldi method;
- - without iterative refinement
- -*/
- -MAT *iter_arnoldi(ip,h_rem,Q,H)
- -ITER *ip;
- -Real *h_rem;
- -MAT *Q, *H;
- -{
- - static VEC *u=VNULL, *r=VNULL, *tmp=VNULL;
- - VEC v; /* auxiliary vector */
- - int i,j;
- - Real h_val, c;
- -
- - if (ip == INULL)
- - error(E_NULL,"iter_arnoldi");
- - if ( ! ip->Ax || ! Q || ! ip->x )
- - error(E_NULL,"iter_arnoldi");
- - if ( ip->k <= 0 )
- - error(E_BOUNDS,"iter_arnoldi");
- - if ( Q->n != ip->x->dim || Q->m != ip->k )
- - error(E_SIZES,"iter_arnoldi");
- -
- - m_zero(Q);
- - H = m_resize(H,ip->k,ip->k);
- - m_zero(H);
- -
- - u = v_resize(u,ip->x->dim);
- - tmp = v_resize(tmp,ip->x->dim);
- - r = v_resize(r,ip->k);
- - MEM_STAT_REG(u,TYPE_VEC);
- - MEM_STAT_REG(tmp,TYPE_VEC);
- - MEM_STAT_REG(r,TYPE_VEC);
- -
- - v.dim = v.max_dim = ip->x->dim;
- -
- - c = v_norm2(ip->x);
- - if ( c <= 0.0)
- - return H;
- - else {
- - v.ve = Q->me[0];
- - sv_mlt(1.0/c,ip->x,&v);
- - }
- -
- - v_zero(r);
- - for ( i = 0; i < ip->k; i++ )
- - {
- - v.ve = Q->me[i];
- - u = (ip->Ax)(ip->A_par,&v,u);
- - v_zero(tmp);
- - for (j = 0; j <= i; j++) {
- - v.ve = Q->me[j];
- - r->ve[j] = in_prod(&v,u);
- - v_mltadd(tmp,&v,r->ve[j],tmp);
- - }
- - v_sub(u,tmp,u);
- - h_val = v_norm2(u);
- - /* if u == 0 then we have an exact subspace */
- - if ( h_val <= 0.0 )
- - {
- - *h_rem = h_val;
- - return H;
- - }
- - set_col(H,i,r);
- - if ( i == ip->k-1 )
- - {
- - *h_rem = h_val;
- - continue;
- - }
- - /* H->me[i+1][i] = h_val; */
- - m_set_val(H,i+1,i,h_val);
- - v.ve = Q->me[i+1];
- - sv_mlt(1.0/h_val,u,&v);
- - }
- -
- - return H;
- -}
- -
- -
- -
- -/* iter_sparnoldi -- uses arnoldi() with an explicit representation of A */
- -MAT *iter_sparnoldi(A,x0,m,h_rem,Q,H)
- -SPMAT *A;
- -VEC *x0;
- -int m;
- -Real *h_rem;
- -MAT *Q, *H;
- -{
- - ITER *ip;
- -
- - ip = iter_get(0,0);
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *) A;
- - ip->x = x0;
- - ip->k = m;
- - iter_arnoldi_iref(ip,h_rem,Q,H);
- - ip->shared_x = ip->shared_b = TRUE;
- - iter_free(ip); /* release only ITER structure */
- - return H;
- -}
- -
- -
- -/* for testing gmres */
- -static void test_gmres(ip,i,Q,R,givc,givs,h_val)
- -ITER *ip;
- -int i;
- -MAT *Q, *R;
- -VEC *givc, *givs;
- -double h_val;
- -{
- - VEC vt, vt1;
- - static MAT *Q1, *R1;
- - int j;
- -
- - /* test Q*A*Q^T = R */
- -
- - Q = m_resize(Q,i+1,ip->b->dim);
- - Q1 = m_resize(Q1,i+1,ip->b->dim);
- - R1 = m_resize(R1,i+1,i+1);
- - MEM_STAT_REG(Q1,TYPE_MAT);
- - MEM_STAT_REG(R1,TYPE_MAT);
- -
- - vt.dim = vt.max_dim = ip->b->dim;
- - vt1.dim = vt1.max_dim = ip->b->dim;
- - for (j=0; j <= i; j++) {
- - vt.ve = Q->me[j];
- - vt1.ve = Q1->me[j];
- - ip->Ax(ip->A_par,&vt,&vt1);
- - }
- -
- - mmtr_mlt(Q,Q1,R1);
- - R1 = m_resize(R1,i+2,i+1);
- - for (j=0; j < i; j++)
- - R1->me[i+1][j] = 0.0;
- - R1->me[i+1][i] = h_val;
- -
- - for (j = 0; j <= i; j++) {
- - rot_rows(R1,j,j+1,givc->ve[j],givs->ve[j],R1);
- - }
- -
- - R1 = m_resize(R1,i+1,i+1);
- - m_sub(R,R1,R1);
- - /* if (m_norm_inf(R1) > MACHEPS*ip->b->dim) */
- - printf(" %d. ||Q*A*Q^T - H|| = %g [cf. MACHEPS = %g]\n",
- - ip->steps,m_norm_inf(R1),MACHEPS);
- -
- - /* check Q*Q^T = I */
- -
- - Q = m_resize(Q,i+1,ip->b->dim);
- - mmtr_mlt(Q,Q,R1);
- - for (j=0; j <= i; j++)
- - R1->me[j][j] -= 1.0;
- - if (m_norm_inf(R1) > MACHEPS*ip->b->dim)
- - printf(" ! m_norm_inf(Q*Q^T) = %g\n",m_norm_inf(R1));
- -
- -}
- -
- -
- -/* gmres -- generalised minimum residual algorithm of Saad & Schultz
- - SIAM J. Sci. Stat. Comp. v.7, pp.856--869 (1986)
- -*/
- -VEC *iter_gmres(ip)
- -ITER *ip;
- -{
- - static VEC *u=VNULL, *r=VNULL, *rhs = VNULL;
- - static VEC *givs=VNULL, *givc=VNULL, *z = VNULL;
- - static MAT *Q = MNULL, *R = MNULL;
- - VEC *rr, v, v1; /* additional pointers (not real vectors) */
- - int i,j, done;
- - Real nres;
- -/* Real last_h; */
- -
- - if (ip == INULL)
- - error(E_NULL,"iter_gmres");
- - if ( ! ip->Ax || ! ip->b )
- - error(E_NULL,"iter_gmres");
- - if ( ! ip->stop_crit )
- - error(E_NULL,"iter_gmres");
- - if ( ip->k <= 0 )
- - error(E_BOUNDS,"iter_gmres");
- - if (ip->x != VNULL && ip->x->dim != ip->b->dim)
- - error(E_SIZES,"iter_gmres");
- -
- - r = v_resize(r,ip->k+1);
- - u = v_resize(u,ip->b->dim);
- - rhs = v_resize(rhs,ip->k+1);
- - givs = v_resize(givs,ip->k); /* Givens rotations */
- - givc = v_resize(givc,ip->k);
- -
- - MEM_STAT_REG(r,TYPE_VEC);
- - MEM_STAT_REG(u,TYPE_VEC);
- - MEM_STAT_REG(rhs,TYPE_VEC);
- - MEM_STAT_REG(givs,TYPE_VEC);
- - MEM_STAT_REG(givc,TYPE_VEC);
- -
- - R = m_resize(R,ip->k+1,ip->k);
- - Q = m_resize(Q,ip->k,ip->b->dim);
- - MEM_STAT_REG(R,TYPE_MAT);
- - MEM_STAT_REG(Q,TYPE_MAT);
- -
- - if (ip->x == VNULL) { /* ip->x == 0 */
- - ip->x = v_get(ip->b->dim);
- - ip->shared_x = FALSE;
- - }
- -
- - v.dim = v.max_dim = ip->b->dim; /* v and v1 are pointers to rows */
- - v1.dim = v1.max_dim = ip->b->dim; /* of matrix Q */
- -
- - if (ip->Bx != (Fun_Ax)NULL) { /* if precondition is defined */
- - z = v_resize(z,ip->b->dim);
- - MEM_STAT_REG(z,TYPE_VEC);
- - }
- -
- - done = FALSE;
- - for (ip->steps = 0; ip->steps <= ip->limit; ) {
- -
- - /* restart */
- -
- - ip->Ax(ip->A_par,ip->x,u); /* u = A*x */
- - v_sub(ip->b,u,u); /* u = b - A*x */
- - rr = u; /* rr is a pointer only */
- -
- - if (ip->Bx) {
- - (ip->Bx)(ip->B_par,u,z); /* tmp = B*(b-A*x) */
- - rr = z;
- - }
- -
- - nres = v_norm2(rr);
- - if (ip->steps == 0) {
- - if (ip->info) ip->info(ip,nres,VNULL,VNULL);
- - if ( ip->stop_crit(ip,nres,VNULL,VNULL) ) {
- - done = TRUE;
- - break;
- - }
- - }
- - else if (nres <= 0.0) break;
- -
- - v.ve = Q->me[0];
- - sv_mlt(1.0/nres,rr,&v);
- -
- - v_zero(r);
- - v_zero(rhs);
- - rhs->ve[0] = nres;
- -
- - for ( i = 0; i < ip->k; i++ ) {
- - ip->steps++;
- - v.ve = Q->me[i];
- - (ip->Ax)(ip->A_par,&v,u);
- - rr = u;
- - if (ip->Bx) {
- - (ip->Bx)(ip->B_par,u,z);
- - rr = z;
- - }
- -
- - if (i < ip->k - 1) {
- - v1.ve = Q->me[i+1];
- - v_copy(rr,&v1);
- - for (j = 0; j <= i; j++) {
- - v.ve = Q->me[j];
- - r->ve[j] = in_prod(&v,rr);
- - v_mltadd(&v1,&v,-r->ve[j],&v1);
- - }
- -
- - r->ve[i+1] = nres = v_norm2(&v1);
- - if (nres <= 0.0) {
- - warning(WARN_RES_LESS_0,"iter_gmres");
- - break;
- - }
- - sv_mlt(1.0/nres,&v1,&v1);
- - }
- - else { /* i == ip->k - 1 */
- - /* Q->me[ip->k] need not be computed */
- -
- - for (j = 0; j <= i; j++) {
- - v.ve = Q->me[j];
- - r->ve[j] = in_prod(&v,rr);
- - }
- -
- - nres = in_prod(rr,rr) - in_prod(r,r);
- - if (nres <= 0.0) {
- - warning(WARN_RES_LESS_0,"iter_gmres");
- - break;
- - }
- - r->ve[i+1] = sqrt(nres);
- - }
- -
- - /* QR update */
- -
- - /* last_h = r->ve[i+1]; */ /* for test only */
- - for (j = 0; j < i; j++)
- - rot_vec(r,j,j+1,givc->ve[j],givs->ve[j],r);
- - givens(r->ve[i],r->ve[i+1],&givc->ve[i],&givs->ve[i]);
- - rot_vec(r,i,i+1,givc->ve[i],givs->ve[i],r);
- - rot_vec(rhs,i,i+1,givc->ve[i],givs->ve[i],rhs);
- -
- - set_col(R,i,r);
- -
- - nres = fabs((double) rhs->ve[i+1]);
- - if (ip->info) ip->info(ip,nres,VNULL,VNULL);
- - if ( ip->stop_crit(ip,nres,VNULL,VNULL) ||
- - ip->steps >= ip->limit ) {
- - done = TRUE;
- - break;
- - }
- - }
- -
- - /* use ixi submatrix of R */
- -
- - if (nres <= 0.0) {
- - i--;
- - done = TRUE;
- - }
- - if (i == ip->k) i = ip->k - 1;
- -
- - R = m_resize(R,i+1,i+1);
- - rhs = v_resize(rhs,i+1);
- -
- - /* test only */
- - /* test_gmres(ip,i,Q,R,givc,givs,last_h); */
- -
- - Usolve(R,rhs,rhs,0.0); /* solve a system: R*x = rhs */
- -
- - /* new approximation */
- -
- - for (j = 0; j <= i; j++) {
- - v.ve = Q->me[j];
- - v_mltadd(ip->x,&v,rhs->ve[j],ip->x);
- - }
- -
- - if (done) break;
- -
- - /* back to old dimensions */
- -
- - rhs = v_resize(rhs,ip->k+1);
- - R = m_resize(R,ip->k+1,ip->k);
- -
- - }
- -
- - return ip->x;
- -}
- -
- -/* iter_spgmres - a simple interface to iter_gmres */
- -
- -VEC *iter_spgmres(A,B,b,tol,x,k,limit,steps)
- -SPMAT *A, *B;
- -VEC *b, *x;
- -double tol;
- -int *steps,k,limit;
- -{
- - ITER *ip;
- -
- - ip = iter_get(0,0);
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *) A;
- - if (B) {
- - ip->Bx = (Fun_Ax) sp_mv_mlt;
- - ip->B_par = (void *) B;
- - }
- - else {
- - ip->Bx = (Fun_Ax) NULL;
- - ip->B_par = NULL;
- - }
- - ip->k = k;
- - ip->limit = limit;
- - ip->info = (Fun_info) NULL;
- - ip->b = b;
- - ip->eps = tol;
- - ip->x = x;
- - iter_gmres(ip);
- - x = ip->x;
- - if (steps) *steps = ip->steps;
- - ip->shared_x = ip->shared_b = TRUE;
- - iter_free(ip); /* release only ITER structure */
- - return x;
- -}
- -
- -
- -/* for testing mgcr */
- -static void test_mgcr(ip,i,Q,R)
- -ITER *ip;
- -int i;
- -MAT *Q, *R;
- -{
- - VEC vt, vt1;
- - static MAT *R1;
- - static VEC *r, *r1;
- - VEC *rr;
- - int k,j;
- - Real sm;
- -
- -
- - /* check Q*Q^T = I */
- - vt.dim = vt.max_dim = ip->b->dim;
- - vt1.dim = vt1.max_dim = ip->b->dim;
- -
- - Q = m_resize(Q,i+1,ip->b->dim);
- - R1 = m_resize(R1,i+1,i+1);
- - r = v_resize(r,ip->b->dim);
- - r1 = v_resize(r1,ip->b->dim);
- - MEM_STAT_REG(R1,TYPE_MAT);
- - MEM_STAT_REG(r,TYPE_VEC);
- - MEM_STAT_REG(r1,TYPE_VEC);
- -
- - m_zero(R1);
- - for (k=1; k <= i; k++)
- - for (j=1; j <= i; j++) {
- - vt.ve = Q->me[k];
- - vt1.ve = Q->me[j];
- - R1->me[k][j] = in_prod(&vt,&vt1);
- - }
- - for (j=1; j <= i; j++)
- - R1->me[j][j] -= 1.0;
- - if (m_norm_inf(R1) > MACHEPS*ip->b->dim)
- - printf(" ! (mgcr:) m_norm_inf(Q*Q^T) = %g\n",m_norm_inf(R1));
- -
- - /* check (r_i,Ap_j) = 0 for j <= i */
- -
- - ip->Ax(ip->A_par,ip->x,r);
- - v_sub(ip->b,r,r);
- - rr = r;
- - if (ip->Bx) {
- - ip->Bx(ip->B_par,r,r1);
- - rr = r1;
- - }
- -
- - printf(" ||r|| = %g\n",v_norm2(rr));
- - sm = 0.0;
- - for (j = 1; j <= i; j++) {
- - vt.ve = Q->me[j];
- - sm = max(sm,in_prod(&vt,rr));
- - }
- - if (sm >= MACHEPS*ip->b->dim)
- - printf(" ! (mgcr:) max_j (r,Ap_j) = %g\n",sm);
- -
- -}
- -
- -
- -
- -
- -/*
- - iter_mgcr -- modified generalized conjugate residual algorithm;
- - fast version of GCR;
- -*/
- -VEC *iter_mgcr(ip)
- -ITER *ip;
- -{
- - static VEC *As, *beta, *alpha, *z;
- - static MAT *N, *H;
- -
- - VEC *rr, v, s; /* additional pointer and structures */
- - Real nres; /* norm of a residual */
- - Real dd; /* coefficient d_i */
- - int i,j;
- - int done; /* if TRUE then stop the iterative process */
- - int dim; /* dimension of the problem */
- -
- - /* ip cannot be NULL */
- - if (ip == INULL) error(E_NULL,"mgcr");
- - /* Ax, b and stopping criterion must be given */
- - if (! ip->Ax || ! ip->b || ! ip->stop_crit)
- - error(E_NULL,"mgcr");
- - /* at least one direction vector must exist */
- - if ( ip->k <= 0) error(E_BOUNDS,"mgcr");
- - /* if the vector x is given then b and x must have the same dimension */
- - if ( ip->x && ip->x->dim != ip->b->dim)
- - error(E_SIZES,"mgcr");
- -
- - dim = ip->b->dim;
- - As = v_resize(As,dim);
- - alpha = v_resize(alpha,ip->k);
- - beta = v_resize(beta,ip->k);
- -
- - MEM_STAT_REG(As,TYPE_VEC);
- - MEM_STAT_REG(alpha,TYPE_VEC);
- - MEM_STAT_REG(beta,TYPE_VEC);
- -
- - H = m_resize(H,ip->k,ip->k);
- - N = m_resize(N,ip->k,dim);
- -
- - MEM_STAT_REG(H,TYPE_MAT);
- - MEM_STAT_REG(N,TYPE_MAT);
- -
- - /* if a preconditioner is defined */
- - if (ip->Bx) {
- - z = v_resize(z,dim);
- - MEM_STAT_REG(z,TYPE_VEC);
- - }
- -
- - /* if x is NULL then it is assumed that x has
- - entries with value zero */
- - if ( ! ip->x ) {
- - ip->x = v_get(ip->b->dim);
- - ip->shared_x = FALSE;
- - }
- -
- - /* v and s are additional pointers to rows of N */
- - /* they must have the same dimension as rows of N */
- - v.dim = v.max_dim = s.dim = s.max_dim = dim;
- -
- -
- - done = FALSE;
- - ip->steps = 0;
- - while (TRUE) {
- - (*ip->Ax)(ip->A_par,ip->x,As); /* As = A*x */
- - v_sub(ip->b,As,As); /* As = b - A*x */
- - rr = As; /* rr is an additional pointer */
- -
- - /* if a preconditioner is defined */
- - if (ip->Bx) {
- - (*ip->Bx)(ip->B_par,As,z); /* z = B*(b-A*x) */
- - rr = z;
- - }
- -
- - /* norm of the residual */
- - nres = v_norm2(rr);
- - dd = nres*nres; /* dd = ||r_i||^2 */
- -
- - /* we need to check if the norm of the residual is small enough
- - only when we start the iterative process;
- - otherwise it has been checked yet.
- - Also the member ip->init_res is updated indirectly by
- - ip->stop_crit.
- - */
- - if (ip->steps == 0) { /* information for a user */
- - if (ip->info) (*ip->info)(ip,nres,As,rr);
- - if ( (*ip->stop_crit)(ip,nres,As,rr) ) {
- - /* iterative process is finished */
- - done = TRUE;
- - break;
- - }
- - }
- - else if (nres <= 0.0)
- - break; /* residual is zero -> finish the process */
- -
- - /* save this residual in the first row of N */
- - v.ve = N->me[0];
- - v_copy(rr,&v);
- -
- - for (i = 0; i < ip->k && ip->steps <= ip->limit; i++) {
- - ip->steps++;
- - v.ve = N->me[i]; /* pointer to a row of N (=s_i) */
- - /* note that we must use here &v, not v */
- - (*ip->Ax)(ip->A_par,&v,As);
- - rr = As; /* As = A*s_i */
- - if (ip->Bx) {
- - (*ip->Bx)(ip->B_par,As,z); /* z = B*A*s_i */
- - rr = z;
- - }
- -
- - if (i < ip->k - 1) {
- - s.ve = N->me[i+1]; /* pointer to a row of N (=s_{i+1}) */
- - v_copy(rr,&s); /* s_{i+1} = B*A*s_i */
- - for (j = 0; j <= i-1; j++) {
- - v.ve = N->me[j+1]; /* pointer to a row of N (=s_{j+1}) */
- - beta->ve[j] = in_prod(&v,rr); /* beta_{j,i} */
- - /* s_{i+1} -= beta_{j,i}*s_{j+1} */
- - v_mltadd(&s,&v,- beta->ve[j],&s);
- - }
- -
- - /* beta_{i,i} = ||s_{i+1}||_2 */
- - beta->ve[i] = nres = v_norm2(&s);
- - if ( nres <= 0.0) break; /* if s_{i+1} == 0 then finish */
- - sv_mlt(1.0/nres,&s,&s); /* normalize s_{i+1} */
- -
- - v.ve = N->me[0];
- - alpha->ve[i] = in_prod(&v,&s); /* alpha_i = (s_0 , s_{i+1}) */
- -
- - }
- - else {
- - for (j = 0; j <= i-1; j++) {
- - v.ve = N->me[j+1]; /* pointer to a row of N (=s_{j+1}) */
- - beta->ve[j] = in_prod(&v,rr); /* beta_{j,i} */
- - }
- -
- - nres = in_prod(rr,rr); /* rr = B*A*s_{k-1} */
- - for (j = 0; j <= i-1; j++)
- - nres -= beta->ve[j]*beta->ve[j];
- - if (nres <= 0.0) break; /* s_k is zero */
- - else beta->ve[i] = sqrt(nres); /* beta_{k-1,k-1} */
- -
- - v.ve = N->me[0];
- - alpha->ve[i] = in_prod(&v,rr);
- - for (j = 0; j <= i-1; j++)
- - alpha->ve[i] -= beta->ve[j]*alpha->ve[j];
- - alpha->ve[i] /= beta->ve[i]; /* alpha_{k-1} */
- -
- - }
- -
- - set_col(H,i,beta);
- -
- - dd -= alpha->ve[i]*alpha->ve[i];
- - nres = sqrt(fabs((double) dd));
- - if (dd < 0.0) { /* do restart */
- - if (ip->info) (*ip->info)(ip,-nres,VNULL,VNULL);
- - break;
- - }
- -
- - if (ip->info) (*ip->info)(ip,nres,VNULL,VNULL);
- - if ( ip->stop_crit(ip,nres,VNULL,VNULL) ) {
- - /* stopping criterion is satisfied */
- - done = TRUE;
- - break;
- - }
- -
- - } /* end of for */
- -
- - if (nres <= 0.0) {
- - i--;
- - done = TRUE;
- - }
- - if (i >= ip->k) i = ip->k - 1;
- -
- - /* use (i+1) by (i+1) submatrix of H */
- - H = m_resize(H,i+1,i+1);
- - alpha = v_resize(alpha,i+1);
- - Usolve(H,alpha,alpha,0.0); /* c_i is saved in alpha */
- -
- - for (j = 0; j <= i; j++) {
- - v.ve = N->me[j];
- - v_mltadd(ip->x,&v,alpha->ve[j],ip->x);
- - }
- -
- -
- - if (done) break; /* stop the iterative process */
- - alpha = v_resize(alpha,ip->k);
- - H = m_resize(H,ip->k,ip->k);
- -
- - } /* end of while */
- -
- - return ip->x; /* return the solution */
- -}
- -
- -
- -
- -/* iter_spmgcr - a simple interface to iter_mgcr */
- -/* no preconditioner */
- -VEC *iter_spmgcr(A,B,b,tol,x,k,limit,steps)
- -SPMAT *A, *B;
- -VEC *b, *x;
- -double tol;
- -int *steps,k,limit;
- -{
- - ITER *ip;
- -
- - ip = iter_get(0,0);
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *) A;
- - if (B) {
- - ip->Bx = (Fun_Ax) sp_mv_mlt;
- - ip->B_par = (void *) B;
- - }
- - else {
- - ip->Bx = (Fun_Ax) NULL;
- - ip->B_par = NULL;
- - }
- -
- - ip->k = k;
- - ip->limit = limit;
- - ip->info = (Fun_info) NULL;
- - ip->b = b;
- - ip->eps = tol;
- - ip->x = x;
- - iter_mgcr(ip);
- - x = ip->x;
- - if (steps) *steps = ip->steps;
- - ip->shared_x = ip->shared_b = TRUE;
- - iter_free(ip); /* release only ITER structure */
- - return x;
- -}
- -
- -
- -
- -/*
- - Conjugate gradients method for a normal equation
- - a preconditioner B must be symmetric !!
- -*/
- -VEC *iter_cgne(ip)
- -ITER *ip;
- -{
- - static VEC *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
- - Real alpha, beta, inner, old_inner, nres;
- - VEC *rr1; /* pointer only */
- -
- - if (ip == INULL)
- - error(E_NULL,"iter_cgne");
- - if (!ip->Ax || ! ip->ATx || !ip->b)
- - error(E_NULL,"iter_cgne");
- - if ( ip->x == ip->b )
- - error(E_INSITU,"iter_cgne");
- - if (!ip->stop_crit)
- - error(E_NULL,"iter_cgne");
- -
- - if ( ip->eps <= 0.0 )
- - ip->eps = MACHEPS;
- -
- - r = v_resize(r,ip->b->dim);
- - p = v_resize(p,ip->b->dim);
- - q = v_resize(q,ip->b->dim);
- -
- - MEM_STAT_REG(r,TYPE_VEC);
- - MEM_STAT_REG(p,TYPE_VEC);
- - MEM_STAT_REG(q,TYPE_VEC);
- -
- - z = v_resize(z,ip->b->dim);
- - MEM_STAT_REG(z,TYPE_VEC);
- -
- - if (ip->x) {
- - if (ip->x->dim != ip->b->dim)
- - error(E_SIZES,"iter_cgne");
- - ip->Ax(ip->A_par,ip->x,p); /* p = A*x */
- - v_sub(ip->b,p,z); /* z = b - A*x */
- - }
- - else { /* ip->x == 0 */
- - ip->x = v_get(ip->b->dim);
- - ip->shared_x = FALSE;
- - v_copy(ip->b,z);
- - }
- - rr1 = z;
- - if (ip->Bx) {
- - (ip->Bx)(ip->B_par,rr1,p);
- - rr1 = p;
- - }
- - (ip->ATx)(ip->AT_par,rr1,r); /* r = A^T*B*(b-A*x) */
- -
- -
- - old_inner = 0.0;
- - for ( ip->steps = 0; ip->steps <= ip->limit; ip->steps++ )
- - {
- - rr1 = r;
- - if ( ip->Bx ) {
- - (ip->Bx)(ip->B_par,r,z); /* rr = B*r */
- - rr1 = z;
- - }
- -
- - inner = in_prod(r,rr1);
- - nres = sqrt(fabs(inner));
- - if (ip->info) ip->info(ip,nres,r,rr1);
- - if ( ip->stop_crit(ip,nres,r,rr1) ) break;
- -
- - if ( ip->steps ) /* if ( ip->steps > 0 ) ... */
- - {
- - beta = inner/old_inner;
- - p = v_mltadd(rr1,p,beta,p);
- - }
- - else /* if ( ip->steps == 0 ) ... */
- - {
- - beta = 0.0;
- - p = v_copy(rr1,p);
- - old_inner = 0.0;
- - }
- - (ip->Ax)(ip->A_par,p,q); /* q = A*p */
- - if (ip->Bx) {
- - (ip->Bx)(ip->B_par,q,z);
- - (ip->ATx)(ip->AT_par,z,q);
- - rr1 = q; /* q = A^T*B*A*p */
- - }
- - else {
- - (ip->ATx)(ip->AT_par,q,z); /* z = A^T*A*p */
- - rr1 = z;
- - }
- -
- - alpha = inner/in_prod(rr1,p);
- - v_mltadd(ip->x,p,alpha,ip->x);
- - v_mltadd(r,rr1,-alpha,r);
- - old_inner = inner;
- - }
- -
- - return ip->x;
- -}
- -
- -/* iter_spcgne -- a simple interface to iter_cgne() which
- - uses sparse matrix data structures
- - -- assumes that B contains an actual preconditioner (or NULL)
- - use always as follows:
- - x = iter_spcgne(A,B,b,eps,x,limit,steps);
- - or
- - x = iter_spcgne(A,B,b,eps,VNULL,limit,steps);
- - In the second case the solution vector is created.
- -*/
- -VEC *iter_spcgne(A,B,b,eps,x,limit,steps)
- -SPMAT *A, *B;
- -VEC *b, *x;
- -double eps;
- -int *steps, limit;
- -{
- - ITER *ip;
- -
- - ip = iter_get(0,0);
- - ip->Ax = (Fun_Ax) sp_mv_mlt;
- - ip->A_par = (void *)A;
- - ip->ATx = (Fun_Ax) sp_vm_mlt;
- - ip->AT_par = (void *)A;
- - if (B) {
- - ip->Bx = (Fun_Ax) sp_mv_mlt;
- - ip->B_par = (void *)B;
- - }
- - else {
- - ip->Bx = (Fun_Ax) NULL;
- - ip->B_par = NULL;
- - }
- - ip->info = (Fun_info) NULL;
- - ip->b = b;
- - ip->eps = eps;
- - ip->limit = limit;
- - ip->x = x;
- - iter_cgne(ip);
- - x = ip->x;
- - if (steps) *steps = ip->steps;
- - ip->shared_x = ip->shared_b = TRUE;
- - iter_free(ip); /* release only ITER structure */
- - return x;
- -}
- -
- -
- -
- -
- //GO.SYSIN DD iternsym.c
-