Annotation of OpenXM/src/hgm/fisher-bingham/src/ko-io.c, Revision 1.1
1.1 ! takayama 1: /*********************************************
! 2: file name: io.c
! 3: Description: Functions for input and output
! 4: the data.
! 5: **********************************************/
! 6: #include<stdio.h>
! 7: #include<stdlib.h>
! 8: #include<stdarg.h>
! 9: #include <gsl/gsl_matrix.h>
! 10: #include <gsl/gsl_vector.h>
! 11:
! 12: double Two = 0.5;
! 13: int flag_diag = 0;
! 14:
! 15: /*
! 16: int input_from_file(char *fname, double ***Xp, double **Yp)
! 17: Reading the values of dim,X, and Y from a file.
! 18: */
! 19: int
! 20: input_from_file(char *fname, double ***Xp, double **Yp)
! 21: {
! 22: FILE *fp;
! 23: int dim, M, i, j;
! 24: double **X, *Xy, *Y;
! 25:
! 26: fp = fopen(fname,"r");
! 27: if (fp == NULL) {
! 28: fprintf(stderr,"%s:File open error.\n", fname);
! 29: return(-1);
! 30: }
! 31:
! 32: fscanf(fp,"%d", &dim);
! 33: #ifdef _SPHERE
! 34: M = dim+1;
! 35: #else
! 36: M = dim;
! 37: #endif
! 38: if (M>0){
! 39: X = (double **)malloc(sizeof(double *) * M); /* memory allocation */
! 40: Xy = (double *)malloc(sizeof(double) * M*M);
! 41: Y = (double *)malloc(sizeof(double) * M);
! 42: for (i=0; i<M; i++) { /* initializing */
! 43: X[i] = Xy + i*M;
! 44: Y[i] = 0.0;
! 45: for (j=0; j<M; j++){
! 46: X[i][j] = 0.0;
! 47: }
! 48: }
! 49: }else{
! 50: return -1;
! 51: }
! 52:
! 53: if(flag_diag){
! 54: for(i=0; i<M; i++)
! 55: fscanf(fp, "%lf", &(X[i][i]));
! 56: for(i=0; i<M; i++) /* X <- (X + X^T)/2 */
! 57: for(j=i+1; j<M; j++)
! 58: X[i][j] = X[j][i] = 0.0;
! 59: }else{
! 60: for(i=0; i<M; i++)
! 61: for(j=0; j<M; j++)
! 62: fscanf(fp, "%lf", &(X[i][j]));
! 63: for(i=0; i<M; i++) /* X <- (X + X^T)/2 */
! 64: for(j=i+1; j<M; j++)
! 65: X[i][j] = X[j][i] = (X[i][j]+X[j][i])/2.0 * Two;
! 66: }
! 67: for(i=0; i<M; i++)
! 68: fscanf(fp, "%lf", &(Y[i]));
! 69:
! 70: fclose(fp);
! 71:
! 72: *Xp = X, *Yp = Y;
! 73: return dim;
! 74: }
! 75:
! 76:
! 77: /*
! 78: Calclate the dimension.
! 79: This function is called by input_from_command only.
! 80: */
! 81: static int caldim(int nv)
! 82: {
! 83: int s=0,n=2;
! 84: if(flag_diag){
! 85: if(nv > 0 && nv%2 == 0 )
! 86: return (nv/2);
! 87: }else{
! 88: while (s<nv)
! 89: if (nv==(s+=n++)) return(n-2);
! 90: }
! 91: return(-1);
! 92: }
! 93:
! 94: /*
! 95: int input_from_command(const int argc,char *argv[],const int optind, double ***Xp, double **Yp)
! 96: Reading the values of X and Y from
! 97: the command-line arguments.
! 98: This function have to be written after
! 99: "getopt loop".
! 100: The variable "getind" is definde the header file
! 101: "getopt.h".
! 102: */
! 103: int
! 104: input_from_command(const int argc,char *argv[],const int optind, double ***Xp, double **Yp)
! 105: {
! 106: double **X, *Xy, *Y;
! 107: int M = caldim(argc - optind);
! 108: int dim;
! 109: int i,j,k;
! 110:
! 111: if (M>0){
! 112: dim=M-1;
! 113: X = (double **)malloc(sizeof(double *) * M); /* memory allocation */
! 114: Xy = (double *)malloc(sizeof(double) * M*M);
! 115: Y = (double *)malloc(sizeof(double) * M);
! 116: for (i=0; i<M; i++) { /* initializing */
! 117: X[i] = Xy + i*M;
! 118: Y[i] = 0.0;
! 119: for (j=0; j<M; j++){
! 120: X[i][j] = 0.0;
! 121: }
! 122: }
! 123: }else{
! 124: return -1;
! 125: }
! 126:
! 127: k = optind;
! 128: if(flag_diag){
! 129: for (i=0; i<M; i++){ /* reading X[][] from argv[] */
! 130: sscanf(argv[k++],"%lf",&(X[i][i]));
! 131: if (k>= argc) break;
! 132: }
! 133: for (i=0; i<M; i++)
! 134: for(j=i+1; j<M; j++)
! 135: X[i][j] = X[j][i] = 0.0;
! 136: }else{
! 137: for (i=0; i<M; i++) { /* reading X[][] from argv[] */
! 138: for (j=i; j<M; j++) {
! 139: sscanf(argv[k],"%lf",&(X[i][j]));
! 140: if (j>i) X[i][j] *= Two;
! 141: X[j][i] = X[i][j];
! 142: k++;
! 143: if (k >= argc) break;
! 144: }
! 145: if (k>= argc) break;
! 146: }
! 147: }
! 148: for (i=0; i<M; i++) { /* reading Y[] from argv[]*/
! 149: sscanf(argv[k],"%lf",&(Y[i]));
! 150: k++;
! 151: if (k >= argc) break;
! 152: }
! 153: *Xp = X; *Yp = Y;
! 154: return dim;
! 155: }
! 156:
! 157: /*
! 158: void free_xy(double **X, double Y)
! 159: Description:
! 160: Free the memory allocated by the function
! 161: input_from_file or input_from_command.
! 162: */
! 163: void
! 164: free_xy(double **X, double *Y)
! 165: {
! 166: free(*X); free(X);
! 167: free(Y);
! 168: }
! 169:
! 170: /*
! 171: void print_vec(int n, int length, char *name, double v[], ...);
! 172: Description:
! 173: This function is for printing vectors which
! 174: has the same length.
! 175: The variable "n" means the number of vectors
! 176: and "length" means the length of vectors.
! 177: example:
! 178: print_vec(4,10,
! 179: "foo:\n", foo
! 180: "bar:\n", bar,
! 181: "baz:\n", baz,
! 182: "qux:\n", qux
! 183: );
! 184: */
! 185: void
! 186: print_vec(int n, int length, ...)
! 187: {
! 188: va_list ap;
! 189: char *s;
! 190: double *v;
! 191: int i,j;
! 192:
! 193: va_start(ap, length);
! 194: for(i= 0; i<n; i++){
! 195: s = va_arg(ap, char*);
! 196: printf("%s",s);
! 197: v = va_arg(ap, double*);
! 198: for(j=0; j<length;j++)
! 199: printf("%10lg ", v[j]);
! 200: printf("\n");
! 201: }
! 202: va_end(ap);
! 203: }
! 204:
! 205: /*
! 206: void print_gsl_vector(gsl_vector *v, char *str)
! 207: Description:
! 208: This function is for printing gsl_vector with
! 209: name.
! 210: */
! 211: void
! 212: print_gsl_vector(gsl_vector *v, char *str)
! 213: {
! 214: int length = v->size;
! 215: int i,j;
! 216:
! 217: printf("%s\n\t", str);
! 218: for(i=0; i<length; i++){
! 219: printf("%10.6f ", gsl_vector_get(v, i));
! 220: }
! 221: putchar('\n');
! 222: }
! 223:
! 224:
! 225:
! 226: /*
! 227: void print_gsl_matrix(gsl_matrix *m, char *str)
! 228: Description:
! 229: This function is for printing gsl_matrix with
! 230: name.
! 231: */
! 232: void
! 233: print_gsl_matrix(gsl_matrix *m, char *str)
! 234: {
! 235: int nrow = m->size1, ncol = m->size2;
! 236: int i,j;
! 237:
! 238: printf("%s\n", str);
! 239: for(i=0; i<nrow; i++){
! 240: putchar('\t');
! 241: for(j=0; j<ncol; j++){
! 242: printf("%10.6f ", gsl_matrix_get(m, i, j));
! 243: }
! 244: putchar('\n');
! 245: }
! 246: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>