home *** CD-ROM | disk | FTP | other *** search
- initialises x so that x->ve[i] == i */
- VEC *v_count(x)
- VEC *x;
- {
- int i;
-
- if ( ! x )
- error(E_NULL,"v_count");
-
- for ( i = 0; i < x->d 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 error(E_SIZES,"m_sub");
- if ( out==(MAT *)NULL || out->m != mat1->m || out->n != mat1->n )
- out = m_resize(out,mat1->m,mat1->n);
- m = mat1->m; n = mat1->n;
- for ( i=0; i<m; i++ )
- {
- __sub__(mat1->me[i],mat2->me[i],out->me[i],(int)n);
- /**************************************************
- for ( j=0; j<n; j++ )
- out->me[i][j] = mat1->me[i][j]-mat2->me[i][j];
- **************************************************/
- }
-
- return (out);
- }
-
- /* m_mlt -- matrix-matrix multiplication */
- MAT *m_mlt(A,B,OUT)
- MAT *A,*B,*OUT;
- {
- u_int i, /* j, */ k, m, n, p;
- Real **A_v, **B_v /*, *B_row, *OUT_row, sum, tmp */;
-
- if ( A==(MAT *)NULL || B==(MAT *)NULL )
- error(E_NULL,"m_mlt");
- if ( A->n != B->m )
- error(E_SIZES,"m_mlt");
- if ( A == OUT || B == OUT )
- error(E_INSITU,"m_mlt");
- m = A->m; n = A->n; p = B->n;
- A_v = A->me; B_v = B->me;
-
- if ( OUT==(MAT *)NULL || OUT->m != A->m || OUT->n != B->n )
- OUT = m_resize(OUT,A->m,B->n);
-
- /****************************************************************
- for ( i=0; i<m; i++ )
- for ( j=0; j<p; j++ )
- {
- sum = 0.0;
- for ( k=0; k<n; k++ )
- sum += A_v[i][k]*B_v[k][j];
- OUT->me[i][j] = sum;
- }
- ****************************************************************/
- m_zero(OUT);
- for ( i=0; i<m; i++ )
- for ( k=0; k<n; k++ )
- {
- if ( A_v[i][k] != 0.0 )
- __mltadd__(OUT->me[i],B_v[k],A_v[i][k],(int)p);
- /**************************************************
- B_row = B_v[k]; OUT_row = OUT->me[i];
- for ( j=0; j<p; j++ )
- (*OUT_row++) += tmp*(*B_row++);
- **************************************************/
- }
-
- return OUT;
- }
-
- /* mmtr_mlt -- matrix-matrix transposed multiplication
- -- A.B^T is returned, and stored in OUT */
- MAT *mmtr_mlt(A,B,OUT)
- MAT *A, *B, *OUT;
- {
- int i, j, limit;
- /* Real *A_row, *B_row, sum; */
-
- if ( ! A || ! B )
- error(E_NULL,"mmtr_mlt");
- if ( A == OUT || B == OUT )
- error(E_INSITU,"mmtr_mlt");
- if ( A->n != B->n )
- error(E_SIZES,"mmtr_mlt");
- if ( ! OUT || OUT->m != A->m || OUT->n != B->m )
- OUT = m_resize(OUT,A->m,B->m);
-
- limit = A->n;
- for ( i = 0; i < A->m; i++ )
- for ( j = 0; j < B->m; j++ )
- {
- OUT->me[i][j] = __ip__(A->me[i],B->me[j],(int)limit);
- /**************************************************
- sum = 0.0;
- A_row = A->me[i];
- B_row = B->me[j];
- for ( k = 0; k < limit; k++ )
- sum += (*A_row++)*(*B_row++);
- OUT->me[i][j] = sum;
- **************************************************/
- }
-
- return OUT;
- }
-
- /* mtrm_mlt -- matrix transposed-matrix multiplication
- -- A^T.B is returned, result stored in OUT */
- MAT *mtrm_mlt(A,B,OUT)
- MAT *A, *B, *OUT;
- {
- int i, k, limit;
- /* Real *B_row, *OUT_row, multiplier; */
-
- if ( ! A || ! B )
- error(E_NULL,"mmtr_mlt");
- if ( A == OUT || B == OUT )
- error(E_INSITU,"mtrm_mlt");
- if ( A->m != B->m )
- error(E_SIZES,"mmtr_mlt");
- if ( ! OUT || OUT->m != A->n || OUT->n != B->n )
- OUT = m_resize(OUT,A->n,B->n);
-
- limit = B->n;
- m_zero(OUT);
- for ( k = 0; k < A->m; k++ )
- for ( i = 0; i < A->n; i++ )
- {
- if ( A->me[k][i] != 0.0 )
- __mltadd__(OUT->me[i],B->me[k],A->me[k][i],(int)limit);
- /**************************************************
- multiplier = A->me[k][i];
- OUT_row = OUT->me[i];
- B_row = B->me[k];
- for ( j = 0; j < limit; j++ )
- *(OUT_row++) += multiplier*(*B_row++);
- **************************************************/
- }
-
- return OUT;
- }
-
- /* mv_mlt -- matrix-vector multiplication
- -- Note: b is treated as a column vector */
- VEC *mv_mlt(A,b,out)
- MAT *A;
- VEC *b,*out;
- {
- u_int i, m, n;
- Real **A_v, *b_v /*, *A_row */;
- /* register Real sum; */
-
- if ( A==(MAT *)NULL || b==(VEC *)NULL )
- error(E_NULL,"mv_mlt");
- if ( A->n != b->dim )
- error(E_SIZES,"mv_mlt");
- if ( b == out )
- error(E_INSITU,"mv_mlt");
- if ( out == (VEC *)NULL || out->dim != A->m )
- out = v_resize(out,A->m);
-
- m = A->m; n = A->n;
- A_v = A->me; b_v = b->ve;
- for ( i=0; i<m; i++ )
- {
- /* for ( j=0; j<n; j++ )
- sum += A_v[i][j]*b_v[j]; */
- out->ve[i] = __ip__(A_v[i],b_v,(int)n);
- /**************************************************
- A_row = A_v[i]; b_v = b->ve;
- for ( j=0; j<n; j++ )
- sum += (*A_row++)*(*b_v++);
- out->ve[i] = sum;
- **************************************************/
- }
-
- return out;
- }
-
- /* sm_mlt -- scalar-matrix multiply -- may be in-situ */
- MAT *sm_mlt(scalar,matrix,out)
- double scalar;
- MAT *matrix,*out;
- {
- u_int m,n,i;
-
- if ( matrix==(MAT *)NULL )
- error(E_NULL,"sm_mlt");
- if ( out==(MAT *)NULL || out->m != matrix->m || out->n != matrix->n )
- out = m_resize(out,matrix->m,matrix->n);
- m = matrix->m; n = matrix->n;
- for ( i=0; i<m; i++ )
- __smlt__(matrix->me[i],(double)scalar,out->me[i],(int)n);
- /**************************************************
- for ( j=0; j<n; j++ )
- out->me[i][j] = scalar*matrix->me[i][j];
- **************************************************/
- return (out);
- }
-
- /* vm_mlt -- vector-matrix multiplication
- -- Note: b is treated as a row vector */
- VEC *vm_mlt(A,b,out)
- MAT *A;
- VEC *b,*out;
- {
- u_int j,m,n;
- /* Real sum,**A_v,*b_v; */
-
- if ( A==(MAT *)NULL || b==(VEC *)NULL )
- error(E_NULL,"vm_mlt");
- if ( A->m != b->dim )
- error(E_SIZES,"vm_mlt");
- if ( b == out )
- error(E_INSITU,"vm_mlt");
- if ( out == (VEC *)NULL || out->dim != A->n )
- out = v_resize(out,A->n);
-
- m = A->m; n = A->n;
-
- v_zero(out);
- for ( j = 0; j < m; j++ )
- if ( b->ve[j] != 0.0 )
- __mltadd__(out->ve,A->me[j],b->ve[j],(int)n);
- /**************************************************
- A_v = A->me; b_v = b->ve;
- for ( j=0; j<n; j++ )
- {
- sum = 0.0;
- for ( i=0; i<m; i++ )
- sum += b_v[i]*A_v[i][j];
- out->ve[j] = sum;
- }
- **************************************************/
-
- return out;
- }
-
- /* m_transp -- transpose matrix */
- MAT *m_transp(in,out)
- MAT *in, *out;
- {
- int i, j;
- int in_situ;
- Real tmp;
-
- if ( in == (MAT *)NULL )
- error(E_NULL,"m_transp");
- if ( in == out && in->n != in->m )
- error(E_INSITU2,"m_transp");
- in_situ = ( in == out );
- if ( out == (MAT *)NULL || out->m != in->n || out->n != in->m )
- out = m_resize(out,in->n,in->m);
-
- if ( ! in_situ )
- for ( i = 0; i < in->m; i++ )
- for ( j = 0; j < in->n; j++ )
- out->me[j][i] = in->me[i][j];
- else
- for ( i = 1; i < in->m; i++ )
- for ( j = 0; j < i; j++ )
- { tmp = in->me[i][j];
- in->me[i][j] = in->me[j][i];
- in->me[j][i] = tmp;
- }
-
- return out;
- }
-
- /* swap_rows -- swaps rows i and j of matrix A upto column lim */
- MAT *swap_rows(A,i,j,lo,hi)
- MAT *A;
- int i, j, lo, hi;
- {
- int k;
- Real **A_me, tmp;
-
- if ( ! A )
- error(E_NULL,"swap_rows");
- if ( i < 0 || j < 0 || i >= A->m || j >= A->m )
- error(E_SIZES,"swap_rows");
- lo = max(0,lo);
- hi = min(hi,A->n-1);
- A_me = A->me;
-
- for ( k = lo; k <= hi; k++ )
- {
- tmp = A_me[k][i];
- A_me[k][i] = A_me[k][j];
- A_me[k][j] = tmp;
- }
- return A;
- }
-
- /* swap_cols -- swap columns i and j of matrix A upto row lim */
- MAT *swap_cols(A,i,j,lo,hi)
- MAT *A;
- int i, j, lo, hi;
- {
- int k;
- Real **A_me, tmp;
-
- if ( ! A )
- error(E_NULL,"swap_cols");
- if ( i < 0 || j < 0 || i >= A->n || j >= A->n )
- error(E_SIZES,"swap_cols");
- lo = max(0,lo);
- hi = min(hi,A->m-1);
- A_me = A->me;
-
- for ( k = lo; k <= hi; k++ )
- {
- tmp = A_me[i][k];
- A_me[i][k] = A_me[j][k];
- A_me[j][k] = tmp;
- }
- return A;
- }
-
- /* ms_mltadd -- matrix-scalar multiply and add
- -- may be in situ
- -- returns out == A1 + s*A2 */
- MAT *ms_mltadd(A1,A2,s,out)
- MAT *A1, *A2, *out;
- double s;
- {
- /* register Real *A1_e, *A2_e, *out_e; */
- /* register int j; */
- int i, m, n;
-
- if ( ! A1 || ! A2 )
- error(E_NULL,"ms_mltadd");
- if ( A1->m != A2->m || A1->n != A2->n )
- error(E_SIZES,"ms_mltadd");
-
- if ( s == 0.0 )
- return m_copy(A1,out);
- if ( s == 1.0 )
- return m_add(A1,A2,out);
-
- tracecatch(out = m_copy(A1,out),"ms_mltadd");
-
- m = A1->m; n = A1->n;
- for ( i = 0; i < m; i++ )
- {
- __mltadd__(out->me[i],A2->me[i],s,(int)n);
- /**************************************************
- A1_e = A1->me[i];
- A2_e = A2->me[i];
- out_e = out->me[i];
- for ( j = 0; j < n; j++ )
- out_e[j] = A1_e[j] + s*A2_e[j];
- **************************************************/
- }
-
- return out;
- }
-
- /* mv_mltadd -- matrix-vector multiply and add
- -- may not be in situ
- -- returns out == v1 + alpha*A*v2 */
- VEC *mv_mltadd(v1,v2,A,alpha,out)
- VEC *v1, *v2, *out;
- MAT *A;
- double alpha;
- {
- /* register int j; */
- int i, m, n;
- Real *v2_ve, *out_ve;
-
- if ( ! v1 || ! v2 || ! A )
- error(E_NULL,"mv_mltadd");
- if ( out == v2 )
- error(E_INSITU,"mv_mltadd");
- if ( v1->dim != A->m || v2->dim != A-> n )
- error(E_SIZES,"mv_mltadd");
-
- tracecatch(out = v_copy(v1,out),"mv_mltadd");
-
- v2_ve = v2->ve; out_ve = out->ve;
- m = A->m; n = A->n;
-
- if ( alpha == 0.0 )
- return out;
-
- for ( i = 0; i < m; i++ )
- {
- out_ve[i] += alpha*__ip__(A->me[i],v2_ve,(int)n);
- /**************************************************
- A_e = A->me[i];
- sum = 0.0;
- for ( j = 0; j < n; j++ )
- sum += A_e[j]*v2_ve[j];
- out_ve[i] = v1->ve[i] + alpha*sum;
- **************************************************/
- }
-
- return out;
- }
-
- /* vm_mltadd -- vector-matrix multiply and add
- -- may not be in situ
- -- returns out' == v1' + v2'*A */
- VEC *vm_mltadd(v1,v2,A,alpha,out)
- VEC *v1, *v2, *out;
- MAT *A;
- double alpha;
- {
- int /* i, */ j, m, n;
- Real tmp, /* *A_e, */ *out_ve;
-
- if ( ! v1 || ! v2 || ! A )
- error(E_NULL,"vm_mltadd");
- if ( v2 == out )
- error(E_INSITU,"vm_mltadd");
- if ( v1->dim != A->n || A->m != v2->dim )
- error(E_SIZES,"vm_mltadd");
-
- tracecatch(out = v_copy(v1,out),"vm_mltadd");
-
- out_ve = out->ve; m = A->m; n = A->n;
- for ( j = 0; j < m; j++ )
- {
- tmp = v2->ve[j]*alpha;
- if ( tmp != 0.0 )
- __mltadd__(out_ve,A->me[j],tmp,(int)n);
- /**************************************************
- A_e = A->me[j];
- for ( i = 0; i < n; i++ )
- out_ve[i] += A_e[i]*tmp;
- **************************************************/
- }
-
- return out;
- }
-
-