diff -Nrc HLib-1.3orig/Library/clusterbasis.c HLib-1.3/Library/clusterbasis.c *** HLib-1.3orig/Library/clusterbasis.c 2004-12-12 17:42:31.000000000 +0100 --- HLib-1.3/Library/clusterbasis.c 2006-03-10 11:38:11.603150064 +0100 *************** *** 318,349 **** } double ! getentry_clusterbasis(pclusterbasis cb, int i) { pclusterbasis *son = cb->son; int sons = cb->sons; double **T = cb->T; double *V = cb->V; double z; ! int j, f, s; assert(i >= 0 && i < cb->n); if(sons > 0) { ! s = f = 0; for(j=0; jn; if(i < s) break; } assert(j < sons); ! addeval_lapack(son[j]->kt, cb->kt, T[j], cb->yt, son[j]->yt); ! z = getentry_clusterbasis(son[j], i-f); } else ! z = ddot_(&cb->kt, V+i, &cb->n, cb->yt, eins_); return z; } --- 318,352 ---- } double ! getentry_clusterbasis(pcclusterbasis cb, const double *yt, int i) { pclusterbasis *son = cb->son; int sons = cb->sons; double **T = cb->T; double *V = cb->V; + double *yt1; double z; ! int j, off, s; assert(i >= 0 && i < cb->n); if(sons > 0) { ! s = off = 0; for(j=0; jn; if(i < s) break; } assert(j < sons); ! yt1 = allocate_vector(son[j]->k); ! eval_lapack(son[j]->k, cb->k, T[j], yt, yt1); ! z = getentry_clusterbasis(son[j], yt1, i-off); ! freemem(yt1); } else ! z = ddot_(&cb->k, V+i, &cb->n, yt, eins_); return z; } diff -Nrc HLib-1.3orig/Library/clusterbasis.h HLib-1.3/Library/clusterbasis.h *** HLib-1.3orig/Library/clusterbasis.h 2004-12-12 17:42:31.000000000 +0100 --- HLib-1.3/Library/clusterbasis.h 2006-03-10 11:38:49.797343664 +0100 *************** *** 116,126 **** void backward_clusterbasis(pclusterbasis cb, double *y); ! /* Evaluate z = (V yt)_i for this cluster the fast way, i.e., by using ! using the transfer matrices T to transform yt into the basis of the ! leaf containing i and then applying V */ double ! getentry_clusterbasis(pclusterbasis cb, int i); /* Check whether this cluster basis is orthogonal, i.e., if V^T V = I holds tol: Tolerance, a basis is considered orthonal if |V^T V - I|_F > tol */ --- 116,128 ---- void backward_clusterbasis(pclusterbasis cb, double *y); ! /* Evaluate z = (V yt)_i by using the transfer matrices T to transform ! yt into the basis of the leaf containing i and then applying V ! cb: Cluster basis ! yt: Input vector ! i: Component */ double ! getentry_clusterbasis(pcclusterbasis cb, const double *yt, int i); /* Check whether this cluster basis is orthogonal, i.e., if V^T V = I holds tol: Tolerance, a basis is considered orthonal if |V^T V - I|_F > tol */ diff -Nrc HLib-1.3orig/Library/uniformmatrix.c HLib-1.3/Library/uniformmatrix.c *** HLib-1.3orig/Library/uniformmatrix.c 2004-12-12 17:42:33.000000000 +0100 --- HLib-1.3/Library/uniformmatrix.c 2006-03-10 11:39:23.777177944 +0100 *************** *** 843,869 **** { pclusterbasis bcol = u->col; pclusterbasis brow = u->row; double z; int i, j; assert(row >= 0 && row < brow->n); assert(col >= 0 && col < bcol->n); ! clear_vector(brow->kt, brow->xt); ! ! for(i=0; ikt; i++) { ! for(j=0; jkt; j++) ! bcol->yt[j] = (i == j ? 1.0 : 0.0); ! ! z = getentry_clusterbasis(bcol, col); ! daxpy_(&brow->kt, &z, u->S + i*bcol->kt, eins_, brow->xt, eins_); } - dcopy_(&brow->kt, brow->xt, eins_, brow->yt, eins_); - - z = getentry_clusterbasis(brow, row); - return z; } --- 843,876 ---- { pclusterbasis bcol = u->col; pclusterbasis brow = u->row; + double *yr, *yc; double z; int i, j; + z = 0.0; + assert(row >= 0 && row < brow->n); assert(col >= 0 && col < bcol->n); ! if(u->S) { ! yr = allocate_vector(brow->k); ! yc = allocate_vector(bcol->k); ! ! for(i=0; ik; i++) { ! for(j=0; jk; j++) ! yc[j] = (i == j ? 1.0 : 0.0); ! ! z = getentry_clusterbasis(bcol, yc, col); ! ! daxpy_(&brow->k, &z, u->S + i*bcol->k, eins_, yr, eins_); ! } ! ! z = getentry_clusterbasis(brow, yr, row); ! freemem(yc); ! freemem(yr); } return z; }