The R Project SVN – R – Diff – /trunk/src/main/array.c



#include "duplicate.h"
 
#include <complex.h>
#include "Rcomplex.h"	/* toC99 */
 
 
 
 
 
 
 
 
 
 
 
/* "GetRowNames" and "GetColNames" are utility routines which
 * locate and return the row names and column names from the
 * dimnames attribute of a matrix.  They are useful because
 * old versions of R used pair-based lists for dimnames
 * whereas recent versions use vector based lists.

    double one = 1.0, zero = 0.0;
    int ione = 1;
 
    if (ncy == 1) /* matrix-vector or dot product */
	F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,
			&nrx, y, &ione, &zero, z, &ione FCONE);
    else if (nrx == 1) /* vector-matrix */
	/* Instead of xY, compute (xY)^T == (Y^T)(x^T)
	   The result is a vector, so transposing its content is no-op */
	F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,
			&nry, x, &ione, &zero, z, &ione FCONE);
    else /* matrix-matrix or outer product */
	F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one, x, 
			&nrx, y, &nry, &zero, z, &nrx FCONE FCONE);
}
 
static void internal_cmatprod(Rcomplex *x, int nrx, int ncx,
                              Rcomplex *y, int nry, int ncy, Rcomplex *z)
{

    char *transa = "N", *transb = "N";
    Rcomplex one, zero;
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
 
    F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
                    x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);
#endif
}
 
static void symcrossprod(double *x, int nr, int nc, double *z)
{

    }
 
    char *trans = "T", *uplo = "U";
    double one = 1.0, zero = 0.0;
 
    F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc
		    FCONE FCONE);
    for (int i = 1; i < nc; i++)
	for (int j = 0; j < i; j++) z[i + NC *j] = z[j + NC * i];
}
 
static void crossprod(double *x, int nrx, int ncx,

    double one = 1.0, zero = 0.0;
    int ione = 1;
 
    if (ncy == 1) /* matrix-vector or dot product */
	F77_CALL(dgemv)(transT, &nrx, &ncx, &one, x,
			&nrx, y, &ione, &zero, z, &ione FCONE);
    else if (ncx == 1) /* vector-matrix */
	/* Instead of (x^T)Y, compute ((x^T)Y)^T == (Y^T)x
	   The result is a vector, so transposing its content is no-op */
	F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,
			&nry, x, &ione, &zero, z, &ione FCONE);
    else /* matrix-matrix  or outer product */
	F77_CALL(dgemm)(transT, transN, &ncx, &ncy, &nrx, &one,
		        x, &nrx, y, &nry, &zero, z, &ncx FCONE FCONE);
}
 
static void ccrossprod(Rcomplex *x, int nrx, int ncx,
		       Rcomplex *y, int nry, int ncy, Rcomplex *z)
{

    char *transa = "T", *transb = "N";
    Rcomplex one, zero;
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
 
    F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
                    x, &nrx, y, &nry, &zero, z, &ncx FCONE FCONE);
#endif
}
 
static void symtcrossprod(double *x, int nr, int nc, double *z)
{

    }
 
    char *trans = "N", *uplo = "U";
    double one = 1.0, zero = 0.0;
 
    F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr
		    FCONE FCONE);
    for (int i = 1; i < nr; i++)
	for (int j = 0; j < i; j++) z[i + nr *j] = z[j + nr * i];
}
 
static void tcrossprod(double *x, int nrx, int ncx,

    double one = 1.0, zero = 0.0;
    int ione = 1;
 
    if (nry == 1) /* matrix-vector or dot product */
	F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,
			&nrx, y, &ione, &zero, z, &ione FCONE);
    else if (nrx == 1) /* vector-matrix */
	/* Instead of x(Y^T), compute (x(Y^T))^T == Y(x^T)
	   The result is a vector, so transposing its content is no-op */
	F77_CALL(dgemv)(transN, &nry, &ncy, &one, y,
			&nry, x, &ione, &zero, z, &ione FCONE);
    else /* matrix-matrix or outer product */
	F77_CALL(dgemm)(transN, transT, &nrx, &nry, &ncx, &one,
			x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);
}
 
static void tccrossprod(Rcomplex *x, int nrx, int ncx,
			Rcomplex *y, int nry, int ncy, Rcomplex *z)
{

    char *transa = "N", *transb = "T";
    Rcomplex one, zero;
    one.r = 1.0; one.i = zero.r = zero.i = 0.0;
 
    F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,
                    x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);
#endif
}
 
 
/* "%*%" (op = 0), crossprod (op = 1) or tcrossprod (op = 2) */

       /* copy (part) cols of b to ans */
	for(R_xlen_t j = 0; j < ncb; j++)
	    memcpy(REAL(ans) + j*k, REAL(b) + j*nrb, (size_t)k *sizeof(double));
	double one = 1.0;
	F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",
			&k, &ncb, &one, rr, &nrr, REAL(ans), &k
			FCONE FCONE FCONE FCONE);
    }
    UNPROTECT(nprot);
    return ans;
}
 

The R Project SVN R

(root)/trunk/src/main/array.c – Rev 76540 → 76585

Rev 76540		Rev 76585
Line 33...		Line 33...
33	`#include "duplicate.h"`	33	`#include "duplicate.h"`
34		34
35	`#include <complex.h>`	35	`#include <complex.h>`
36	`#include "Rcomplex.h" /* toC99 */`	36	`#include "Rcomplex.h" /* toC99 */`
37		37
38	`/* FIXME`	-
39	`This calls BLAS routines`	-
40		-
41	`dgemm dgemv zgemm dsyrk dtrsm`	-
42		-
43	`all of which have character args so gfortran (at least) expects a`	-
44	`length passed -- but non-Fortran BLAS may not. At least with cdecl`	-
45	`calling conventions additional arg(s) will not matter.`	-
46	`*/`	-
47		-
48	`/* "GetRowNames" and "GetColNames" are utility routines which`	38	`/* "GetRowNames" and "GetColNames" are utility routines which`
49	`* locate and return the row names and column names from the`	39	`* locate and return the row names and column names from the`
50	`* dimnames attribute of a matrix. They are useful because`	40	`* dimnames attribute of a matrix. They are useful because`
51	`* old versions of R used pair-based lists for dimnames`	41	`* old versions of R used pair-based lists for dimnames`
52	`* whereas recent versions use vector based lists.`	42	`* whereas recent versions use vector based lists.`
Line 809...		Line 799...
809	`double one = 1.0, zero = 0.0;`	799	`double one = 1.0, zero = 0.0;`
810	`int ione = 1;`	800	`int ione = 1;`
811		801
812	`if (ncy == 1) /* matrix-vector or dot product */`	802	`if (ncy == 1) /* matrix-vector or dot product */`
813	`F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,`	803	`F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,`
814	`&nrx, y, &ione, &zero, z, &ione);`	804	`&nrx, y, &ione, &zero, z, &ione FCONE);`
815	`else if (nrx == 1) /* vector-matrix */`	805	`else if (nrx == 1) /* vector-matrix */`
816	`/* Instead of xY, compute (xY)^T == (Y^T)(x^T)`	806	`/* Instead of xY, compute (xY)^T == (Y^T)(x^T)`
817	`The result is a vector, so transposing its content is no-op */`	807	`The result is a vector, so transposing its content is no-op */`
818	`F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,`	808	`F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,`
819	`&nry, x, &ione, &zero, z, &ione);`	809	`&nry, x, &ione, &zero, z, &ione FCONE);`
820	`else /* matrix-matrix or outer product */`	810	`else /* matrix-matrix or outer product */`
821	`F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one,`	811	`F77_CALL(dgemm)(transN, transN, &nrx, &ncy, &ncx, &one, x,`
822	`x, &nrx, y, &nry, &zero, z, &nrx);`	812	`&nrx, y, &nry, &zero, z, &nrx FCONE FCONE);`
823	`}`	813	`}`
824		814
825	`static void internal_cmatprod(Rcomplex *x, int nrx, int ncx,`	815	`static void internal_cmatprod(Rcomplex *x, int nrx, int ncx,`
826	`Rcomplex y, int nry, int ncy, Rcomplex z)`	816	`Rcomplex y, int nry, int ncy, Rcomplex z)`
827	`{`	817	`{`
Line 955...		Line 945...
955	`char transa = "N", transb = "N";`	945	`char transa = "N", transb = "N";`
956	`Rcomplex one, zero;`	946	`Rcomplex one, zero;`
957	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`	947	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`
958		948
959	`F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,`	949	`F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,`
960	`x, &nrx, y, &nry, &zero, z, &nrx);`	950	`x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);`
961	`#endif`	951	`#endif`
962	`}`	952	`}`
963		953
964	`static void symcrossprod(double x, int nr, int nc, double z)`	954	`static void symcrossprod(double x, int nr, int nc, double z)`
965	`{`	955	`{`
Line 992...		Line 982...
992	`}`	982	`}`
993		983
994	`char trans = "T", uplo = "U";`	984	`char trans = "T", uplo = "U";`
995	`double one = 1.0, zero = 0.0;`	985	`double one = 1.0, zero = 0.0;`
996		986
997	`F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);`	987	`F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc`
-		988	`FCONE FCONE);`
998	`for (int i = 1; i < nc; i++)`	989	`for (int i = 1; i < nc; i++)`
999	`for (int j = 0; j < i; j++) z[i + NC j] = z[j + NC i];`	990	`for (int j = 0; j < i; j++) z[i + NC j] = z[j + NC i];`
1000	`}`	991	`}`
1001		992
1002	`static void crossprod(double *x, int nrx, int ncx,`	993	`static void crossprod(double *x, int nrx, int ncx,`
Line 1036...		Line 1027...
1036	`double one = 1.0, zero = 0.0;`	1027	`double one = 1.0, zero = 0.0;`
1037	`int ione = 1;`	1028	`int ione = 1;`
1038		1029
1039	`if (ncy == 1) /* matrix-vector or dot product */`	1030	`if (ncy == 1) /* matrix-vector or dot product */`
1040	`F77_CALL(dgemv)(transT, &nrx, &ncx, &one, x,`	1031	`F77_CALL(dgemv)(transT, &nrx, &ncx, &one, x,`
1041	`&nrx, y, &ione, &zero, z, &ione);`	1032	`&nrx, y, &ione, &zero, z, &ione FCONE);`
1042	`else if (ncx == 1) /* vector-matrix */`	1033	`else if (ncx == 1) /* vector-matrix */`
1043	`/* Instead of (x^T)Y, compute ((x^T)Y)^T == (Y^T)x`	1034	`/* Instead of (x^T)Y, compute ((x^T)Y)^T == (Y^T)x`
1044	`The result is a vector, so transposing its content is no-op */`	1035	`The result is a vector, so transposing its content is no-op */`
1045	`F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,`	1036	`F77_CALL(dgemv)(transT, &nry, &ncy, &one, y,`
1046	`&nry, x, &ione, &zero, z, &ione);`	1037	`&nry, x, &ione, &zero, z, &ione FCONE);`
1047	`else /* matrix-matrix or outer product */`	1038	`else /* matrix-matrix or outer product */`
1048	`F77_CALL(dgemm)(transT, transN, &ncx, &ncy, &nrx, &one,`	1039	`F77_CALL(dgemm)(transT, transN, &ncx, &ncy, &nrx, &one,`
1049	`x, &nrx, y, &nry, &zero, z, &ncx);`	1040	`x, &nrx, y, &nry, &zero, z, &ncx FCONE FCONE);`
1050	`}`	1041	`}`
1051		1042
1052	`static void ccrossprod(Rcomplex *x, int nrx, int ncx,`	1043	`static void ccrossprod(Rcomplex *x, int nrx, int ncx,`
1053	`Rcomplex y, int nry, int ncy, Rcomplex z)`	1044	`Rcomplex y, int nry, int ncy, Rcomplex z)`
1054	`{`	1045	`{`
Line 1090...		Line 1081...
1090	`char transa = "T", transb = "N";`	1081	`char transa = "T", transb = "N";`
1091	`Rcomplex one, zero;`	1082	`Rcomplex one, zero;`
1092	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`	1083	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`
1093		1084
1094	`F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,`	1085	`F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,`
1095	`x, &nrx, y, &nry, &zero, z, &ncx);`	1086	`x, &nrx, y, &nry, &zero, z, &ncx FCONE FCONE);`
1096	`#endif`	1087	`#endif`
1097	`}`	1088	`}`
1098		1089
1099	`static void symtcrossprod(double x, int nr, int nc, double z)`	1090	`static void symtcrossprod(double x, int nr, int nc, double z)`
1100	`{`	1091	`{`
Line 1127...		Line 1118...
1127	`}`	1118	`}`
1128		1119
1129	`char trans = "N", uplo = "U";`	1120	`char trans = "N", uplo = "U";`
1130	`double one = 1.0, zero = 0.0;`	1121	`double one = 1.0, zero = 0.0;`
1131		1122
1132	`F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr);`	1123	`F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr`
-		1124	`FCONE FCONE);`
1133	`for (int i = 1; i < nr; i++)`	1125	`for (int i = 1; i < nr; i++)`
1134	`for (int j = 0; j < i; j++) z[i + nr j] = z[j + nr i];`	1126	`for (int j = 0; j < i; j++) z[i + nr j] = z[j + nr i];`
1135	`}`	1127	`}`
1136		1128
1137	`static void tcrossprod(double *x, int nrx, int ncx,`	1129	`static void tcrossprod(double *x, int nrx, int ncx,`
Line 1169...		Line 1161...
1169	`double one = 1.0, zero = 0.0;`	1161	`double one = 1.0, zero = 0.0;`
1170	`int ione = 1;`	1162	`int ione = 1;`
1171		1163
1172	`if (nry == 1) /* matrix-vector or dot product */`	1164	`if (nry == 1) /* matrix-vector or dot product */`
1173	`F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,`	1165	`F77_CALL(dgemv)(transN, &nrx, &ncx, &one, x,`
1174	`&nrx, y, &ione, &zero, z, &ione);`	1166	`&nrx, y, &ione, &zero, z, &ione FCONE);`
1175	`else if (nrx == 1) /* vector-matrix */`	1167	`else if (nrx == 1) /* vector-matrix */`
1176	`/* Instead of x(Y^T), compute (x(Y^T))^T == Y(x^T)`	1168	`/* Instead of x(Y^T), compute (x(Y^T))^T == Y(x^T)`
1177	`The result is a vector, so transposing its content is no-op */`	1169	`The result is a vector, so transposing its content is no-op */`
1178	`F77_CALL(dgemv)(transN, &nry, &ncy, &one, y,`	1170	`F77_CALL(dgemv)(transN, &nry, &ncy, &one, y,`
1179	`&nry, x, &ione, &zero, z, &ione);`	1171	`&nry, x, &ione, &zero, z, &ione FCONE);`
1180	`else /* matrix-matrix or outer product */`	1172	`else /* matrix-matrix or outer product */`
1181	`F77_CALL(dgemm)(transN, transT, &nrx, &nry, &ncx, &one,`	1173	`F77_CALL(dgemm)(transN, transT, &nrx, &nry, &ncx, &one,`
1182	`x, &nrx, y, &nry, &zero, z, &nrx);`	1174	`x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);`
1183	`}`	1175	`}`
1184		1176
1185	`static void tccrossprod(Rcomplex *x, int nrx, int ncx,`	1177	`static void tccrossprod(Rcomplex *x, int nrx, int ncx,`
1186	`Rcomplex y, int nry, int ncy, Rcomplex z)`	1178	`Rcomplex y, int nry, int ncy, Rcomplex z)`
1187	`{`	1179	`{`
Line 1222...		Line 1214...
1222	`char transa = "N", transb = "T";`	1214	`char transa = "N", transb = "T";`
1223	`Rcomplex one, zero;`	1215	`Rcomplex one, zero;`
1224	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`	1216	`one.r = 1.0; one.i = zero.r = zero.i = 0.0;`
1225		1217
1226	`F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,`	1218	`F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,`
1227	`x, &nrx, y, &nry, &zero, z, &nrx);`	1219	`x, &nrx, y, &nry, &zero, z, &nrx FCONE FCONE);`
1228	`#endif`	1220	`#endif`
1229	`}`	1221	`}`
1230		1222
1231		1223
1232	`/* "%%" (op = 0), crossprod (op = 1) or tcrossprod (op = 2) /`	1224	`/* "%%" (op = 0), crossprod (op = 1) or tcrossprod (op = 2) /`
Line 2291...		Line 2283...
2291	`/* copy (part) cols of b to ans */`	2283	`/* copy (part) cols of b to ans */`
2292	`for(R_xlen_t j = 0; j < ncb; j++)`	2284	`for(R_xlen_t j = 0; j < ncb; j++)`
2293	`memcpy(REAL(ans) + jk, REAL(b) + jnrb, (size_t)k *sizeof(double));`	2285	`memcpy(REAL(ans) + jk, REAL(b) + jnrb, (size_t)k *sizeof(double));`
2294	`double one = 1.0;`	2286	`double one = 1.0;`
2295	`F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",`	2287	`F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",`
2296	`&k, &ncb, &one, rr, &nrr, REAL(ans), &k);`	2288	`&k, &ncb, &one, rr, &nrr, REAL(ans), &k`
-		2289	`FCONE FCONE FCONE FCONE);`
2297	`}`	2290	`}`
2298	`UNPROTECT(nprot);`	2291	`UNPROTECT(nprot);`
2299	`return ans;`	2292	`return ans;`
2300	`}`	2293	`}`
2301		2294