mlfma.c
/* Recompress an ACA approximation using truncated singular values. */
Recompress an ACA approximation using truncated singular values.
Recompress an ACA approximation using truncated singular values .
static int recaca (cplx *u, cplx *v, int m, int n, int k, real tol) {
	int rank, lwork, info, i, j, l;
	real *rwork, *ss;
	cplx *qu, *qv, *rp, *work, *tu, *tv, *us, *vs, alpha = 1.0, beta = 0.0;
	qu = malloc (k * (m + n + 3 * k) * sizeof(cplx));
	qv = qu + k * m;
	rp = qv + k * n;
	us = rp + k * k;
	vs = us + k * k;
	tu = malloc (2 * k * sizeof(cplx));
	tv = tu + k;
	ss = malloc (6 * k * sizeof(real));
	rwork = ss + k;
	lwork = -1;
	GEQRF (&m, &k, qu, &m, tu, qu, &lwork, &info);
	GEQRF (&n, &k, qv, &n, tv, qu + 1, &lwork, &info);
	UNGQR (&m, &k, &k, qu, &m, tu, qu + 2, &lwork, &info);
	UNGQR (&n, &k, &k, qv, &n, tv, qu + 3, &lwork, &info);
	GESVD ("S", "S", &k, &k, rp, &k, ss, us, &k, vs, &k, qu + 4, &lwork, rwork, &info);
	for (i = 0; i < 5; ++i) lwork = MAX(lwork, (int)creal(qu[i]));
	work = malloc(lwork * sizeof(cplx));
	memcpy (qu, u, (m * k) * sizeof(cplx));
	memcpy (qv, v, (n * k) * sizeof(cplx));
	GEQRF (&m, &k, qu, &m, tu, work, &lwork, &info);
	GEQRF (&n, &k, qv, &n, tv, work, &lwork, &info);
	for (i = 0; i < k; ++i) {
		for (j = 0; j < k; ++j) {
			rp[i + j * k] = 0.0;
			for (l = MAX(i,j); l < k; ++l)
				rp[i + j * k] += qu[i + l * m] * conj(qv[j + l * n]);
		}
	}
	GESVD ("S", "S", &k, &k, rp, &k, ss, us, &k, vs, &k, work, &lwork, rwork, &info);
	for (rank = 0; rank < k; ++rank)
		if (fabs(ss[rank] / ss[0]) < tol) break;
	for (i = 0; i < rank; ++i)
		for (j = 0; j < k; ++j)
			vs[i + k * j] *= ss[i];
	UNGQR (&m, &k, &k, qu, &m, tu, work, &lwork, &info);
	UNGQR (&n, &k, &k, qv, &n, tv, work, &lwork, &info);
	GEMM (CblasColMajor, CblasNoTrans, CblasNoTrans, m, rank, k,
			&alpha, qu, m, us, k, &beta, u, m);
	GEMM (CblasColMajor, CblasNoTrans, CblasConjTrans, n, rank, k,
			&alpha, qv, n, vs, k, &beta, v, n);
	free (qu);
	free (tu);
	free (ss);
	free (work);
	return rank;
}
static int recaca ( cplx * u , cplx * v , int m , int n , int k , real tol ) { int rank , lwork , info , i , j , l ; real * rwork , * ss ; cplx * qu , * qv , * rp , * work , * tu , * tv , * us , * vs , alpha = 1.0 , beta = 0.0 ; qu = malloc ( k * ( m + n + 3 * k ) * sizeof ( cplx ) ) ; qv = qu + k * m ; rp = qv + k * n ; us = rp + k * k ; vs = us + k * k ; tu = malloc ( 2 * k * sizeof ( cplx ) ) ; tv = tu + k ; ss = malloc ( 6 * k * sizeof ( real ) ) ; rwork = ss + k ; lwork = -1 ; GEQRF ( & m , & k , qu , & m , tu , qu , & lwork , & info ) ; GEQRF ( & n , & k , qv , & n , tv , qu + 1 , & lwork , & info ) ; UNGQR ( & m , & k , & k , qu , & m , tu , qu + 2 , & lwork , & info ) ; UNGQR ( & n , & k , & k , qv , & n , tv , qu + 3 , & lwork , & info ) ; GESVD ( " " , " " , & k , & k , rp , & k , ss , us , & k , vs , & k , qu + 4 , & lwork , rwork , & info ) ; for ( i = 0 ; i < 5 ; ++ i ) lwork = MAX ( lwork , ( int ) creal ( qu [ i ] ) ) ; work = malloc ( lwork * sizeof ( cplx ) ) ; memcpy ( qu , u , ( m * k ) * sizeof ( cplx ) ) ; memcpy ( qv , v , ( n * k ) * sizeof ( cplx ) ) ; GEQRF ( & m , & k , qu , & m , tu , work , & lwork , & info ) ; GEQRF ( & n , & k , qv , & n , tv , work , & lwork , & info ) ; for ( i = 0 ; i < k ; ++ i ) { for ( j = 0 ; j < k ; ++ j ) { rp [ i + j * k ] = 0.0 ; for ( l = MAX ( i , j ) ; l < k ; ++ l ) rp [ i + j * k ] += qu [ i + l * m ] * conj ( qv [ j + l * n ] ) ; } } GESVD ( " " , " " , & k , & k , rp , & k , ss , us , & k , vs , & k , work , & lwork , rwork , & info ) ; for ( rank = 0 ; rank < k ; ++ rank ) if ( fabs ( ss [ rank ] / ss [ 0 ] ) < tol ) break ; for ( i = 0 ; i < rank ; ++ i ) for ( j = 0 ; j < k ; ++ j ) vs [ i + k * j ] *= ss [ i ] ; UNGQR ( & m , & k , & k , qu , & m , tu , work , & lwork , & info ) ; UNGQR ( & n , & k , & k , qv , & n , tv , work , & lwork , & info ) ; GEMM ( CblasColMajor , CblasNoTrans , CblasNoTrans , m , rank , k , & alpha , qu , m , us , k , & beta , u , m ) ; GEMM ( CblasColMajor , CblasNoTrans , CblasConjTrans , n , rank , k , & alpha , qv , n , vs , k , & beta , v , n ) ; free ( qu ) ; free ( tu ) ; free ( ss ) ; free ( work ) ; return rank ; }
Recompress an ACA approximation using truncated singular values.
Recompress an ACA approximation using truncated singular values .
/* Allocate space for the orthogonal matrices and the outer product. */ /* Allocate space for the reflector coefficients. */ /* Allocate space for the singular values and real work. */ /* Figure out the maximum size of the work array. */ /* Find the maximum workspace size. */ /* Copy the input matrices into the workspaces. */ /* Compute the QR factorizations of the row and column matrices. */ /* Compute the outer product of the triangular matrices. */ /* Compute the SVD of the triangular matrix product. */ /* Find the truncated rank. */ /* Multiply the singular values by the column matrix. */ /* Expand the orthogonal matrices. */ /* Compute the new row and column matrices. */