zlarfb.c

来自「算断裂的」· C语言代码 · 共 736 行 · 第 1/2 页
736 行

                W := C2 */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    zcopy_(m, &C(1,*n-*k+j), &c__1, &WORK(1,j), &c__1);
/* L100: */
		}

/*              W := W * V2 */

		ztrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b1, 
			&V(*n-*k+1,1), ldv, &WORK(1,1), 
			ldwork);
		if (*n > *k) {

/*                 W := W + C1 * V1 */

		    i__1 = *n - *k;
		    zgemm_("No transpose", "No transpose", m, k, &i__1, &c_b1,
			     &C(1,1), ldc, &V(1,1), ldv, &c_b1, &
			    WORK(1,1), ldwork);
		}

/*              W := W * T  or  W * T' */

		ztrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b1, &T(1,1), ldt, &WORK(1,1), ldwork);

/*              C := C - W * V' */

		if (*n > *k) {

/*                 C1 := C1 - W * V1' */

		    i__1 = *n - *k;
		    z__1.r = -1., z__1.i = 0.;
		    zgemm_("No transpose", "Conjugate transpose", m, &i__1, k,
			     &z__1, &WORK(1,1), ldwork, &V(1,1), 
			    ldv, &c_b1, &C(1,1), ldc);
		}

/*              W := W * V2' */

		ztrmm_("Right", "Upper", "Conjugate transpose", "Unit", m, k, 
			&c_b1, &V(*n-*k+1,1), ldv, &WORK(1,1), ldwork);

/*              C2 := C2 - W */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			i__3 = i + (*n - *k + j) * c_dim1;
			i__4 = i + (*n - *k + j) * c_dim1;
			i__5 = i + j * work_dim1;
			z__1.r = C(i,*n-*k+j).r - WORK(i,j).r, z__1.i = C(i,*n-*k+j).i 
				- WORK(i,j).i;
			C(i,*n-*k+j).r = z__1.r, C(i,*n-*k+j).i = z__1.i;
/* L110: */
		    }
/* L120: */
		}
	    }
	}

    } else if (lsame_(storev, "R")) {

	if (lsame_(direct, "F")) {

/*           Let  V =  ( V1  V2 )    (V1: first K columns)   
             where  V1  is unit upper triangular. */

	    if (lsame_(side, "L")) {

/*              Form  H * C  or  H' * C  where  C = ( C1 )   
                                                    ( C2 )   

                W := C' * V'  =  (C1'*V1' + C2'*V2') (stored i
n WORK)   

                W := C1' */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    zcopy_(n, &C(j,1), ldc, &WORK(1,j), &
			    c__1);
		    zlacgv_(n, &WORK(1,j), &c__1);
/* L130: */
		}

/*              W := W * V1' */

		ztrmm_("Right", "Upper", "Conjugate transpose", "Unit", n, k, 
			&c_b1, &V(1,1), ldv, &WORK(1,1), ldwork);
		if (*m > *k) {

/*                 W := W + C2'*V2' */

		    i__1 = *m - *k;
		    zgemm_("Conjugate transpose", "Conjugate transpose", n, k,
			     &i__1, &c_b1, &C(*k+1,1), ldc, &V(1,*k+1), ldv, &c_b1, &WORK(1,1), 
			    ldwork);
		}

/*              W := W * T'  or  W * T */

		ztrmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b1, &T(1,1), ldt, &WORK(1,1), ldwork);

/*              C := C - V' * W' */

		if (*m > *k) {

/*                 C2 := C2 - V2' * W' */

		    i__1 = *m - *k;
		    z__1.r = -1., z__1.i = 0.;
		    zgemm_("Conjugate transpose", "Conjugate transpose", &
			    i__1, n, k, &z__1, &V(1,*k+1), ldv,
			     &WORK(1,1), ldwork, &c_b1, &C(*k+1,1), ldc);
		}

/*              W := W * V1 */

		ztrmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b1, 
			&V(1,1), ldv, &WORK(1,1), ldwork);

/*              C1 := C1 - W' */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    i__2 = *n;
		    for (i = 1; i <= *n; ++i) {
			i__3 = j + i * c_dim1;
			i__4 = j + i * c_dim1;
			d_cnjg(&z__2, &WORK(i,j));
			z__1.r = C(j,i).r - z__2.r, z__1.i = C(j,i).i - 
				z__2.i;
			C(j,i).r = z__1.r, C(j,i).i = z__1.i;
/* L140: */
		    }
/* L150: */
		}

	    } else if (lsame_(side, "R")) {

/*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
   

                W := C * V'  =  (C1*V1' + C2*V2')  (stored in 
WORK)   

                W := C1 */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    zcopy_(m, &C(1,j), &c__1, &WORK(1,j), &c__1);
/* L160: */
		}

/*              W := W * V1' */

		ztrmm_("Right", "Upper", "Conjugate transpose", "Unit", m, k, 
			&c_b1, &V(1,1), ldv, &WORK(1,1), ldwork);
		if (*n > *k) {

/*                 W := W + C2 * V2' */

		    i__1 = *n - *k;
		    zgemm_("No transpose", "Conjugate transpose", m, k, &i__1,
			     &c_b1, &C(1,*k+1), ldc, &V(1,*k+1), ldv, &c_b1, &WORK(1,1), 
			    ldwork);
		}

/*              W := W * T  or  W * T' */

		ztrmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b1, &T(1,1), ldt, &WORK(1,1), ldwork);

/*              C := C - W * V */

		if (*n > *k) {

/*                 C2 := C2 - W * V2 */

		    i__1 = *n - *k;
		    z__1.r = -1., z__1.i = 0.;
		    zgemm_("No transpose", "No transpose", m, &i__1, k, &z__1,
			     &WORK(1,1), ldwork, &V(1,*k+1), ldv, &c_b1, &C(1,*k+1), ldc);
		}

/*              W := W * V1 */

		ztrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b1, 
			&V(1,1), ldv, &WORK(1,1), ldwork);

/*              C1 := C1 - W */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			i__3 = i + j * c_dim1;
			i__4 = i + j * c_dim1;
			i__5 = i + j * work_dim1;
			z__1.r = C(i,j).r - WORK(i,j).r, z__1.i = C(i,j).i 
				- WORK(i,j).i;
			C(i,j).r = z__1.r, C(i,j).i = z__1.i;
/* L170: */
		    }
/* L180: */
		}

	    }

	} else {

/*           Let  V =  ( V1  V2 )    (V2: last K columns)   
             where  V2  is unit lower triangular. */

	    if (lsame_(side, "L")) {

/*              Form  H * C  or  H' * C  where  C = ( C1 )   
                                                    ( C2 )   

                W := C' * V'  =  (C1'*V1' + C2'*V2') (stored i
n WORK)   

                W := C2' */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    zcopy_(n, &C(*m-*k+j,1), ldc, &WORK(1,j), &c__1);
		    zlacgv_(n, &WORK(1,j), &c__1);
/* L190: */
		}

/*              W := W * V2' */

		ztrmm_("Right", "Lower", "Conjugate transpose", "Unit", n, k, 
			&c_b1, &V(1,*m-*k+1), ldv, &WORK(1,1), ldwork);
		if (*m > *k) {

/*                 W := W + C1'*V1' */

		    i__1 = *m - *k;
		    zgemm_("Conjugate transpose", "Conjugate transpose", n, k,
			     &i__1, &c_b1, &C(1,1), ldc, &V(1,1), 
			    ldv, &c_b1, &WORK(1,1), ldwork);
		}

/*              W := W * T'  or  W * T */

		ztrmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b1, &T(1,1), ldt, &WORK(1,1), ldwork);

/*              C := C - V' * W' */

		if (*m > *k) {

/*                 C1 := C1 - V1' * W' */

		    i__1 = *m - *k;
		    z__1.r = -1., z__1.i = 0.;
		    zgemm_("Conjugate transpose", "Conjugate transpose", &
			    i__1, n, k, &z__1, &V(1,1), ldv, &WORK(1,1), ldwork, &c_b1, &C(1,1), ldc);
		}

/*              W := W * V2 */

		ztrmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b1, 
			&V(1,*m-*k+1), ldv, &WORK(1,1), ldwork);

/*              C2 := C2 - W' */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    i__2 = *n;
		    for (i = 1; i <= *n; ++i) {
			i__3 = *m - *k + j + i * c_dim1;
			i__4 = *m - *k + j + i * c_dim1;
			d_cnjg(&z__2, &WORK(i,j));
			z__1.r = C(*m-*k+j,i).r - z__2.r, z__1.i = C(*m-*k+j,i).i - 
				z__2.i;
			C(*m-*k+j,i).r = z__1.r, C(*m-*k+j,i).i = z__1.i;
/* L200: */
		    }
/* L210: */
		}

	    } else if (lsame_(side, "R")) {

/*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
   

                W := C * V'  =  (C1*V1' + C2*V2')  (stored in 
WORK)   

                W := C2 */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    zcopy_(m, &C(1,*n-*k+j), &c__1, &WORK(1,j), &c__1);
/* L220: */
		}

/*              W := W * V2' */

		ztrmm_("Right", "Lower", "Conjugate transpose", "Unit", m, k, 
			&c_b1, &V(1,*n-*k+1), ldv, &WORK(1,1), ldwork);
		if (*n > *k) {

/*                 W := W + C1 * V1' */

		    i__1 = *n - *k;
		    zgemm_("No transpose", "Conjugate transpose", m, k, &i__1,
			     &c_b1, &C(1,1), ldc, &V(1,1), ldv, &
			    c_b1, &WORK(1,1), ldwork);
		}

/*              W := W * T  or  W * T' */

		ztrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b1, &T(1,1), ldt, &WORK(1,1), ldwork);

/*              C := C - W * V */

		if (*n > *k) {

/*                 C1 := C1 - W * V1 */

		    i__1 = *n - *k;
		    z__1.r = -1., z__1.i = 0.;
		    zgemm_("No transpose", "No transpose", m, &i__1, k, &z__1,
			     &WORK(1,1), ldwork, &V(1,1), ldv, &
			    c_b1, &C(1,1), ldc);
		}

/*              W := W * V2 */

		ztrmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b1, 
			&V(1,*n-*k+1), ldv, &WORK(1,1), ldwork);

/*              C1 := C1 - W */

		i__1 = *k;
		for (j = 1; j <= *k; ++j) {
		    i__2 = *m;
		    for (i = 1; i <= *m; ++i) {
			i__3 = i + (*n - *k + j) * c_dim1;
			i__4 = i + (*n - *k + j) * c_dim1;
			i__5 = i + j * work_dim1;
			z__1.r = C(i,*n-*k+j).r - WORK(i,j).r, z__1.i = C(i,*n-*k+j).i 
				- WORK(i,j).i;
			C(i,*n-*k+j).r = z__1.r, C(i,*n-*k+j).i = z__1.i;
/* L230: */
		    }
/* L240: */
		}

	    }

	}
    }

    return 0;

/*     End of ZLARFB */

} /* zlarfb_ */
zlarfb.c - 源码说明

本页面展示了「算断裂的」中的 zlarfb.c 源码文件，采用 C语言编程语言编写，共 736 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。
虫虫下载站收录了大量与断裂相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。
⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?