1 2 /* 3 Module Name: xyt 4 Module Info: 5 6 author: Henry M. Tufo III 7 e-mail: hmt@asci.uchicago.edu 8 contact: 9 +--------------------------------+--------------------------------+ 10 |MCS Division - Building 221 |Department of Computer Science | 11 |Argonne National Laboratory |Ryerson 152 | 12 |9700 S. Cass Avenue |The University of Chicago | 13 |Argonne, IL 60439 |Chicago, IL 60637 | 14 |(630) 252-5354/5986 ph/fx |(773) 702-6019/8487 ph/fx | 15 +--------------------------------+--------------------------------+ 16 17 Last Modification: 3.20.01 18 */ 19 #include <../src/ksp/pc/impls/tfs/tfs.h> 20 21 #define LEFT -1 22 #define RIGHT 1 23 #define BOTH 0 24 25 typedef struct xyt_solver_info { 26 PetscInt n, m, n_global, m_global; 27 PetscInt nnz, max_nnz, msg_buf_sz; 28 PetscInt *nsep, *lnsep, *fo, nfo, *stages; 29 PetscInt *xcol_sz, *xcol_indices; 30 PetscScalar **xcol_vals, *x, *solve_uu, *solve_w; 31 PetscInt *ycol_sz, *ycol_indices; 32 PetscScalar **ycol_vals, *y; 33 PetscInt nsolves; 34 PetscScalar tot_solve_time; 35 } xyt_info; 36 37 typedef struct matvec_info { 38 PetscInt n, m, n_global, m_global; 39 PetscInt *local2global; 40 PCTFS_gs_ADT PCTFS_gs_handle; 41 PetscErrorCode (*matvec)(struct matvec_info *, PetscScalar *, PetscScalar *); 42 void *grid_data; 43 } mv_info; 44 45 struct xyt_CDT { 46 PetscInt id; 47 PetscInt ns; 48 PetscInt level; 49 xyt_info *info; 50 mv_info *mvi; 51 }; 52 53 static PetscInt n_xyt = 0; 54 static PetscInt n_xyt_handles = 0; 55 56 /* prototypes */ 57 static PetscErrorCode do_xyt_solve(xyt_ADT xyt_handle, PetscScalar *rhs); 58 static PetscErrorCode check_handle(xyt_ADT xyt_handle); 59 static PetscErrorCode det_separators(xyt_ADT xyt_handle); 60 static PetscErrorCode do_matvec(mv_info *A, PetscScalar *v, PetscScalar *u); 61 static PetscErrorCode xyt_generate(xyt_ADT xyt_handle); 62 static PetscErrorCode do_xyt_factor(xyt_ADT xyt_handle); 63 static mv_info *set_mvi(PetscInt *local2global, PetscInt n, PetscInt m, PetscErrorCode (*matvec)(mv_info *, PetscScalar *, PetscScalar *), void *grid_data); 64 65 xyt_ADT XYT_new(void) { 66 xyt_ADT xyt_handle; 67 68 /* rolling count on n_xyt ... pot. problem here */ 69 n_xyt_handles++; 70 xyt_handle = (xyt_ADT)malloc(sizeof(struct xyt_CDT)); 71 xyt_handle->id = ++n_xyt; 72 xyt_handle->info = NULL; 73 xyt_handle->mvi = NULL; 74 75 return (xyt_handle); 76 } 77 78 PetscErrorCode XYT_factor(xyt_ADT xyt_handle, /* prev. allocated xyt handle */ 79 PetscInt *local2global, /* global column mapping */ 80 PetscInt n, /* local num rows */ 81 PetscInt m, /* local num cols */ 82 PetscErrorCode (*matvec)(void *, PetscScalar *, PetscScalar *), /* b_loc=A_local.x_loc */ 83 void *grid_data) /* grid data for matvec */ 84 { 85 PCTFS_comm_init(); 86 check_handle(xyt_handle); 87 88 /* only 2^k for now and all nodes participating */ 89 PetscCheck((1 << (xyt_handle->level = PCTFS_i_log2_num_nodes)) == PCTFS_num_nodes, PETSC_COMM_SELF, PETSC_ERR_PLIB, "only 2^k for now and MPI_COMM_WORLD!!! %d != %d", 1 << PCTFS_i_log2_num_nodes, PCTFS_num_nodes); 90 91 /* space for X info */ 92 xyt_handle->info = (xyt_info *)malloc(sizeof(xyt_info)); 93 94 /* set up matvec handles */ 95 xyt_handle->mvi = set_mvi(local2global, n, m, (PetscErrorCode(*)(mv_info *, PetscScalar *, PetscScalar *))matvec, grid_data); 96 97 /* matrix is assumed to be of full rank */ 98 /* LATER we can reset to indicate rank def. */ 99 xyt_handle->ns = 0; 100 101 /* determine separators and generate firing order - NB xyt info set here */ 102 det_separators(xyt_handle); 103 104 return (do_xyt_factor(xyt_handle)); 105 } 106 107 PetscErrorCode XYT_solve(xyt_ADT xyt_handle, PetscScalar *x, PetscScalar *b) { 108 PCTFS_comm_init(); 109 check_handle(xyt_handle); 110 111 /* need to copy b into x? */ 112 if (b) PCTFS_rvec_copy(x, b, xyt_handle->mvi->n); 113 return do_xyt_solve(xyt_handle, x); 114 } 115 116 PetscErrorCode XYT_free(xyt_ADT xyt_handle) { 117 PCTFS_comm_init(); 118 check_handle(xyt_handle); 119 n_xyt_handles--; 120 121 free(xyt_handle->info->nsep); 122 free(xyt_handle->info->lnsep); 123 free(xyt_handle->info->fo); 124 free(xyt_handle->info->stages); 125 free(xyt_handle->info->solve_uu); 126 free(xyt_handle->info->solve_w); 127 free(xyt_handle->info->x); 128 free(xyt_handle->info->xcol_vals); 129 free(xyt_handle->info->xcol_sz); 130 free(xyt_handle->info->xcol_indices); 131 free(xyt_handle->info->y); 132 free(xyt_handle->info->ycol_vals); 133 free(xyt_handle->info->ycol_sz); 134 free(xyt_handle->info->ycol_indices); 135 free(xyt_handle->info); 136 free(xyt_handle->mvi->local2global); 137 PCTFS_gs_free(xyt_handle->mvi->PCTFS_gs_handle); 138 free(xyt_handle->mvi); 139 free(xyt_handle); 140 141 /* if the check fails we nuke */ 142 /* if NULL pointer passed to free we nuke */ 143 /* if the calls to free fail that's not my problem */ 144 return (0); 145 } 146 147 /* This function is currently not used */ 148 PetscErrorCode XYT_stats(xyt_ADT xyt_handle) { 149 PetscInt op[] = {NON_UNIFORM, GL_MIN, GL_MAX, GL_ADD, GL_MIN, GL_MAX, GL_ADD, GL_MIN, GL_MAX, GL_ADD}; 150 PetscInt fop[] = {NON_UNIFORM, GL_MIN, GL_MAX, GL_ADD}; 151 PetscInt vals[9], work[9]; 152 PetscScalar fvals[3], fwork[3]; 153 154 PetscFunctionBegin; 155 PetscCall(PCTFS_comm_init()); 156 PetscCall(check_handle(xyt_handle)); 157 158 /* if factorization not done there are no stats */ 159 if (!xyt_handle->info || !xyt_handle->mvi) { 160 if (!PCTFS_my_id) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "XYT_stats() :: no stats available!\n")); 161 PetscFunctionReturn(0); 162 } 163 164 vals[0] = vals[1] = vals[2] = xyt_handle->info->nnz; 165 vals[3] = vals[4] = vals[5] = xyt_handle->mvi->n; 166 vals[6] = vals[7] = vals[8] = xyt_handle->info->msg_buf_sz; 167 PetscCall(PCTFS_giop(vals, work, PETSC_STATIC_ARRAY_LENGTH(op) - 1, op)); 168 169 fvals[0] = fvals[1] = fvals[2] = xyt_handle->info->tot_solve_time / xyt_handle->info->nsolves++; 170 PetscCall(PCTFS_grop(fvals, fwork, PETSC_STATIC_ARRAY_LENGTH(fop) - 1, fop)); 171 172 if (!PCTFS_my_id) { 173 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: min xyt_nnz=%" PetscInt_FMT "\n", PCTFS_my_id, vals[0])); 174 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: max xyt_nnz=%" PetscInt_FMT "\n", PCTFS_my_id, vals[1])); 175 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: avg xyt_nnz=%g\n", PCTFS_my_id, (double)(1.0 * vals[2] / PCTFS_num_nodes))); 176 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: tot xyt_nnz=%" PetscInt_FMT "\n", PCTFS_my_id, vals[2])); 177 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: xyt C(2d) =%g\n", PCTFS_my_id, (double)(vals[2] / (PetscPowReal(1.0 * vals[5], 1.5))))); 178 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: xyt C(3d) =%g\n", PCTFS_my_id, (double)(vals[2] / (PetscPowReal(1.0 * vals[5], 1.6667))))); 179 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: min xyt_n =%" PetscInt_FMT "\n", PCTFS_my_id, vals[3])); 180 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: max xyt_n =%" PetscInt_FMT "\n", PCTFS_my_id, vals[4])); 181 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: avg xyt_n =%g\n", PCTFS_my_id, (double)(1.0 * vals[5] / PCTFS_num_nodes))); 182 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: tot xyt_n =%" PetscInt_FMT "\n", PCTFS_my_id, vals[5])); 183 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: min xyt_buf=%" PetscInt_FMT "\n", PCTFS_my_id, vals[6])); 184 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: max xyt_buf=%" PetscInt_FMT "\n", PCTFS_my_id, vals[7])); 185 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: avg xyt_buf=%g\n", PCTFS_my_id, (double)(1.0 * vals[8] / PCTFS_num_nodes))); 186 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: min xyt_slv=%g\n", PCTFS_my_id, (double)PetscRealPart(fvals[0]))); 187 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: max xyt_slv=%g\n", PCTFS_my_id, (double)PetscRealPart(fvals[1]))); 188 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d :: avg xyt_slv=%g\n", PCTFS_my_id, (double)PetscRealPart(fvals[2] / PCTFS_num_nodes))); 189 } 190 PetscFunctionReturn(0); 191 } 192 193 /* 194 195 Description: get A_local, local portion of global coarse matrix which 196 is a row dist. nxm matrix w/ n<m. 197 o my_ml holds address of ML struct associated w/A_local and coarse grid 198 o local2global holds global number of column i (i=0,...,m-1) 199 o local2global holds global number of row i (i=0,...,n-1) 200 o mylocmatvec performs A_local . vec_local (note that gs is performed using 201 PCTFS_gs_init/gop). 202 203 mylocmatvec = my_ml->Amat[grid_tag].matvec->external; 204 mylocmatvec (void :: void *data, double *in, double *out) 205 */ 206 static PetscErrorCode do_xyt_factor(xyt_ADT xyt_handle) { 207 return xyt_generate(xyt_handle); 208 } 209 210 static PetscErrorCode xyt_generate(xyt_ADT xyt_handle) { 211 PetscInt i, j, k, idx; 212 PetscInt dim, col; 213 PetscScalar *u, *uu, *v, *z, *w, alpha, alpha_w; 214 PetscInt *segs; 215 PetscInt op[] = {GL_ADD, 0}; 216 PetscInt off, len; 217 PetscScalar *x_ptr, *y_ptr; 218 PetscInt *iptr, flag; 219 PetscInt start = 0, end, work; 220 PetscInt op2[] = {GL_MIN, 0}; 221 PCTFS_gs_ADT PCTFS_gs_handle; 222 PetscInt *nsep, *lnsep, *fo; 223 PetscInt a_n = xyt_handle->mvi->n; 224 PetscInt a_m = xyt_handle->mvi->m; 225 PetscInt *a_local2global = xyt_handle->mvi->local2global; 226 PetscInt level; 227 PetscInt n, m; 228 PetscInt *xcol_sz, *xcol_indices, *stages; 229 PetscScalar **xcol_vals, *x; 230 PetscInt *ycol_sz, *ycol_indices; 231 PetscScalar **ycol_vals, *y; 232 PetscInt n_global; 233 PetscInt xt_nnz = 0, xt_max_nnz = 0; 234 PetscInt yt_nnz = 0, yt_max_nnz = 0; 235 PetscBLASInt i1 = 1, dlen; 236 PetscScalar dm1 = -1.0; 237 238 n = xyt_handle->mvi->n; 239 nsep = xyt_handle->info->nsep; 240 lnsep = xyt_handle->info->lnsep; 241 fo = xyt_handle->info->fo; 242 end = lnsep[0]; 243 level = xyt_handle->level; 244 PCTFS_gs_handle = xyt_handle->mvi->PCTFS_gs_handle; 245 246 /* is there a null space? */ 247 /* LATER add in ability to detect null space by checking alpha */ 248 for (i = 0, j = 0; i <= level; i++) j += nsep[i]; 249 250 m = j - xyt_handle->ns; 251 if (m != j) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "xyt_generate() :: null space exists %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, j, xyt_handle->ns)); 252 253 PetscCall(PetscInfo(0, "xyt_generate() :: X(%" PetscInt_FMT ",%" PetscInt_FMT ")\n", n, m)); 254 255 /* get and initialize storage for x local */ 256 /* note that x local is nxm and stored by columns */ 257 xcol_sz = (PetscInt *)malloc(m * sizeof(PetscInt)); 258 xcol_indices = (PetscInt *)malloc((2 * m + 1) * sizeof(PetscInt)); 259 xcol_vals = (PetscScalar **)malloc(m * sizeof(PetscScalar *)); 260 for (i = j = 0; i < m; i++, j += 2) { 261 xcol_indices[j] = xcol_indices[j + 1] = xcol_sz[i] = -1; 262 xcol_vals[i] = NULL; 263 } 264 xcol_indices[j] = -1; 265 266 /* get and initialize storage for y local */ 267 /* note that y local is nxm and stored by columns */ 268 ycol_sz = (PetscInt *)malloc(m * sizeof(PetscInt)); 269 ycol_indices = (PetscInt *)malloc((2 * m + 1) * sizeof(PetscInt)); 270 ycol_vals = (PetscScalar **)malloc(m * sizeof(PetscScalar *)); 271 for (i = j = 0; i < m; i++, j += 2) { 272 ycol_indices[j] = ycol_indices[j + 1] = ycol_sz[i] = -1; 273 ycol_vals[i] = NULL; 274 } 275 ycol_indices[j] = -1; 276 277 /* size of separators for each sub-hc working from bottom of tree to top */ 278 /* this looks like nsep[]=segments */ 279 stages = (PetscInt *)malloc((level + 1) * sizeof(PetscInt)); 280 segs = (PetscInt *)malloc((level + 1) * sizeof(PetscInt)); 281 PCTFS_ivec_zero(stages, level + 1); 282 PCTFS_ivec_copy(segs, nsep, level + 1); 283 for (i = 0; i < level; i++) segs[i + 1] += segs[i]; 284 stages[0] = segs[0]; 285 286 /* temporary vectors */ 287 u = (PetscScalar *)malloc(n * sizeof(PetscScalar)); 288 z = (PetscScalar *)malloc(n * sizeof(PetscScalar)); 289 v = (PetscScalar *)malloc(a_m * sizeof(PetscScalar)); 290 uu = (PetscScalar *)malloc(m * sizeof(PetscScalar)); 291 w = (PetscScalar *)malloc(m * sizeof(PetscScalar)); 292 293 /* extra nnz due to replication of vertices across separators */ 294 for (i = 1, j = 0; i <= level; i++) j += nsep[i]; 295 296 /* storage for sparse x values */ 297 n_global = xyt_handle->info->n_global; 298 xt_max_nnz = yt_max_nnz = (PetscInt)(2.5 * PetscPowReal(1.0 * n_global, 1.6667) + j * n / 2) / PCTFS_num_nodes; 299 x = (PetscScalar *)malloc(xt_max_nnz * sizeof(PetscScalar)); 300 y = (PetscScalar *)malloc(yt_max_nnz * sizeof(PetscScalar)); 301 302 /* LATER - can embed next sep to fire in gs */ 303 /* time to make the donuts - generate X factor */ 304 for (dim = i = j = 0; i < m; i++) { 305 /* time to move to the next level? */ 306 while (i == segs[dim]) { 307 PetscCheck(dim != level, PETSC_COMM_SELF, PETSC_ERR_PLIB, "dim about to exceed level"); 308 stages[dim++] = i; 309 end += lnsep[dim]; 310 } 311 stages[dim] = i; 312 313 /* which column are we firing? */ 314 /* i.e. set v_l */ 315 /* use new seps and do global min across hc to determine which one to fire */ 316 (start < end) ? (col = fo[start]) : (col = INT_MAX); 317 PCTFS_giop_hc(&col, &work, 1, op2, dim); 318 319 /* shouldn't need this */ 320 if (col == INT_MAX) { 321 PetscCall(PetscInfo(0, "hey ... col==INT_MAX??\n")); 322 continue; 323 } 324 325 /* do I own it? I should */ 326 PCTFS_rvec_zero(v, a_m); 327 if (col == fo[start]) { 328 start++; 329 idx = PCTFS_ivec_linear_search(col, a_local2global, a_n); 330 if (idx != -1) { 331 v[idx] = 1.0; 332 j++; 333 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "NOT FOUND!"); 334 } else { 335 idx = PCTFS_ivec_linear_search(col, a_local2global, a_m); 336 if (idx != -1) v[idx] = 1.0; 337 } 338 339 /* perform u = A.v_l */ 340 PCTFS_rvec_zero(u, n); 341 do_matvec(xyt_handle->mvi, v, u); 342 343 /* uu = X^T.u_l (local portion) */ 344 /* technically only need to zero out first i entries */ 345 /* later turn this into an XYT_solve call ? */ 346 PCTFS_rvec_zero(uu, m); 347 y_ptr = y; 348 iptr = ycol_indices; 349 for (k = 0; k < i; k++) { 350 off = *iptr++; 351 len = *iptr++; 352 PetscCall(PetscBLASIntCast(len, &dlen)); 353 PetscCallBLAS("BLASdot", uu[k] = BLASdot_(&dlen, u + off, &i1, y_ptr, &i1)); 354 y_ptr += len; 355 } 356 357 /* uu = X^T.u_l (comm portion) */ 358 PetscCall(PCTFS_ssgl_radd(uu, w, dim, stages)); 359 360 /* z = X.uu */ 361 PCTFS_rvec_zero(z, n); 362 x_ptr = x; 363 iptr = xcol_indices; 364 for (k = 0; k < i; k++) { 365 off = *iptr++; 366 len = *iptr++; 367 PetscCall(PetscBLASIntCast(len, &dlen)); 368 PetscCallBLAS("BLASaxpy", BLASaxpy_(&dlen, &uu[k], x_ptr, &i1, z + off, &i1)); 369 x_ptr += len; 370 } 371 372 /* compute v_l = v_l - z */ 373 PCTFS_rvec_zero(v + a_n, a_m - a_n); 374 PetscCall(PetscBLASIntCast(n, &dlen)); 375 PetscCallBLAS("BLASaxpy", BLASaxpy_(&dlen, &dm1, z, &i1, v, &i1)); 376 377 /* compute u_l = A.v_l */ 378 if (a_n != a_m) PCTFS_gs_gop_hc(PCTFS_gs_handle, v, "+\0", dim); 379 PCTFS_rvec_zero(u, n); 380 do_matvec(xyt_handle->mvi, v, u); 381 382 /* compute sqrt(alpha) = sqrt(u_l^T.u_l) - local portion */ 383 PetscCall(PetscBLASIntCast(n, &dlen)); 384 PetscCallBLAS("BLASdot", alpha = BLASdot_(&dlen, u, &i1, u, &i1)); 385 /* compute sqrt(alpha) = sqrt(u_l^T.u_l) - comm portion */ 386 PCTFS_grop_hc(&alpha, &alpha_w, 1, op, dim); 387 388 alpha = (PetscScalar)PetscSqrtReal((PetscReal)alpha); 389 390 /* check for small alpha */ 391 /* LATER use this to detect and determine null space */ 392 PetscCheck(PetscAbsScalar(alpha) >= 1.0e-14, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bad alpha! %g", (double)PetscAbsScalar(alpha)); 393 394 /* compute v_l = v_l/sqrt(alpha) */ 395 PCTFS_rvec_scale(v, 1.0 / alpha, n); 396 PCTFS_rvec_scale(u, 1.0 / alpha, n); 397 398 /* add newly generated column, v_l, to X */ 399 flag = 1; 400 off = len = 0; 401 for (k = 0; k < n; k++) { 402 if (v[k] != 0.0) { 403 len = k; 404 if (flag) { 405 off = k; 406 flag = 0; 407 } 408 } 409 } 410 411 len -= (off - 1); 412 413 if (len > 0) { 414 if ((xt_nnz + len) > xt_max_nnz) { 415 PetscCall(PetscInfo(0, "increasing space for X by 2x!\n")); 416 xt_max_nnz *= 2; 417 x_ptr = (PetscScalar *)malloc(xt_max_nnz * sizeof(PetscScalar)); 418 PCTFS_rvec_copy(x_ptr, x, xt_nnz); 419 free(x); 420 x = x_ptr; 421 x_ptr += xt_nnz; 422 } 423 xt_nnz += len; 424 PCTFS_rvec_copy(x_ptr, v + off, len); 425 426 xcol_indices[2 * i] = off; 427 xcol_sz[i] = xcol_indices[2 * i + 1] = len; 428 xcol_vals[i] = x_ptr; 429 } else { 430 xcol_indices[2 * i] = 0; 431 xcol_sz[i] = xcol_indices[2 * i + 1] = 0; 432 xcol_vals[i] = x_ptr; 433 } 434 435 /* add newly generated column, u_l, to Y */ 436 flag = 1; 437 off = len = 0; 438 for (k = 0; k < n; k++) { 439 if (u[k] != 0.0) { 440 len = k; 441 if (flag) { 442 off = k; 443 flag = 0; 444 } 445 } 446 } 447 448 len -= (off - 1); 449 450 if (len > 0) { 451 if ((yt_nnz + len) > yt_max_nnz) { 452 PetscCall(PetscInfo(0, "increasing space for Y by 2x!\n")); 453 yt_max_nnz *= 2; 454 y_ptr = (PetscScalar *)malloc(yt_max_nnz * sizeof(PetscScalar)); 455 PCTFS_rvec_copy(y_ptr, y, yt_nnz); 456 free(y); 457 y = y_ptr; 458 y_ptr += yt_nnz; 459 } 460 yt_nnz += len; 461 PCTFS_rvec_copy(y_ptr, u + off, len); 462 463 ycol_indices[2 * i] = off; 464 ycol_sz[i] = ycol_indices[2 * i + 1] = len; 465 ycol_vals[i] = y_ptr; 466 } else { 467 ycol_indices[2 * i] = 0; 468 ycol_sz[i] = ycol_indices[2 * i + 1] = 0; 469 ycol_vals[i] = y_ptr; 470 } 471 } 472 473 /* close off stages for execution phase */ 474 while (dim != level) { 475 stages[dim++] = i; 476 PetscCall(PetscInfo(0, "disconnected!!! dim(%" PetscInt_FMT ")!=level(%" PetscInt_FMT ")\n", dim, level)); 477 } 478 stages[dim] = i; 479 480 xyt_handle->info->n = xyt_handle->mvi->n; 481 xyt_handle->info->m = m; 482 xyt_handle->info->nnz = xt_nnz + yt_nnz; 483 xyt_handle->info->max_nnz = xt_max_nnz + yt_max_nnz; 484 xyt_handle->info->msg_buf_sz = stages[level] - stages[0]; 485 xyt_handle->info->solve_uu = (PetscScalar *)malloc(m * sizeof(PetscScalar)); 486 xyt_handle->info->solve_w = (PetscScalar *)malloc(m * sizeof(PetscScalar)); 487 xyt_handle->info->x = x; 488 xyt_handle->info->xcol_vals = xcol_vals; 489 xyt_handle->info->xcol_sz = xcol_sz; 490 xyt_handle->info->xcol_indices = xcol_indices; 491 xyt_handle->info->stages = stages; 492 xyt_handle->info->y = y; 493 xyt_handle->info->ycol_vals = ycol_vals; 494 xyt_handle->info->ycol_sz = ycol_sz; 495 xyt_handle->info->ycol_indices = ycol_indices; 496 497 free(segs); 498 free(u); 499 free(v); 500 free(uu); 501 free(z); 502 free(w); 503 504 return (0); 505 } 506 507 static PetscErrorCode do_xyt_solve(xyt_ADT xyt_handle, PetscScalar *uc) { 508 PetscInt off, len, *iptr; 509 PetscInt level = xyt_handle->level; 510 PetscInt n = xyt_handle->info->n; 511 PetscInt m = xyt_handle->info->m; 512 PetscInt *stages = xyt_handle->info->stages; 513 PetscInt *xcol_indices = xyt_handle->info->xcol_indices; 514 PetscInt *ycol_indices = xyt_handle->info->ycol_indices; 515 PetscScalar *x_ptr, *y_ptr, *uu_ptr; 516 PetscScalar *solve_uu = xyt_handle->info->solve_uu; 517 PetscScalar *solve_w = xyt_handle->info->solve_w; 518 PetscScalar *x = xyt_handle->info->x; 519 PetscScalar *y = xyt_handle->info->y; 520 PetscBLASInt i1 = 1, dlen; 521 522 PetscFunctionBegin; 523 uu_ptr = solve_uu; 524 PCTFS_rvec_zero(uu_ptr, m); 525 526 /* x = X.Y^T.b */ 527 /* uu = Y^T.b */ 528 for (y_ptr = y, iptr = ycol_indices; *iptr != -1; y_ptr += len) { 529 off = *iptr++; 530 len = *iptr++; 531 PetscCall(PetscBLASIntCast(len, &dlen)); 532 PetscCallBLAS("BLASdot", *uu_ptr++ = BLASdot_(&dlen, uc + off, &i1, y_ptr, &i1)); 533 } 534 535 /* comunication of beta */ 536 uu_ptr = solve_uu; 537 if (level) PetscCall(PCTFS_ssgl_radd(uu_ptr, solve_w, level, stages)); 538 PCTFS_rvec_zero(uc, n); 539 540 /* x = X.uu */ 541 for (x_ptr = x, iptr = xcol_indices; *iptr != -1; x_ptr += len) { 542 off = *iptr++; 543 len = *iptr++; 544 PetscCall(PetscBLASIntCast(len, &dlen)); 545 PetscCallBLAS("BLASaxpy", BLASaxpy_(&dlen, uu_ptr++, x_ptr, &i1, uc + off, &i1)); 546 } 547 PetscFunctionReturn(0); 548 } 549 550 static PetscErrorCode check_handle(xyt_ADT xyt_handle) { 551 PetscInt vals[2], work[2], op[] = {NON_UNIFORM, GL_MIN, GL_MAX}; 552 553 PetscFunctionBegin; 554 PetscCheck(xyt_handle, PETSC_COMM_SELF, PETSC_ERR_PLIB, "check_handle() :: bad handle :: NULL %p", xyt_handle); 555 556 vals[0] = vals[1] = xyt_handle->id; 557 PCTFS_giop(vals, work, PETSC_STATIC_ARRAY_LENGTH(op) - 1, op); 558 PetscCheck(!(vals[0] != vals[1]) && !(xyt_handle->id <= 0), PETSC_COMM_SELF, PETSC_ERR_PLIB, "check_handle() :: bad handle :: id mismatch min/max %" PetscInt_FMT "/%" PetscInt_FMT " %" PetscInt_FMT, vals[0], vals[1], xyt_handle->id); 559 PetscFunctionReturn(0); 560 } 561 562 static PetscErrorCode det_separators(xyt_ADT xyt_handle) { 563 PetscInt i, ct, id; 564 PetscInt mask, edge, *iptr; 565 PetscInt *dir, *used; 566 PetscInt sum[4], w[4]; 567 PetscScalar rsum[4], rw[4]; 568 PetscInt op[] = {GL_ADD, 0}; 569 PetscScalar *lhs, *rhs; 570 PetscInt *nsep, *lnsep, *fo, nfo = 0; 571 PCTFS_gs_ADT PCTFS_gs_handle = xyt_handle->mvi->PCTFS_gs_handle; 572 PetscInt *local2global = xyt_handle->mvi->local2global; 573 PetscInt n = xyt_handle->mvi->n; 574 PetscInt m = xyt_handle->mvi->m; 575 PetscInt level = xyt_handle->level; 576 PetscInt shared = 0; 577 578 PetscFunctionBegin; 579 dir = (PetscInt *)malloc(sizeof(PetscInt) * (level + 1)); 580 nsep = (PetscInt *)malloc(sizeof(PetscInt) * (level + 1)); 581 lnsep = (PetscInt *)malloc(sizeof(PetscInt) * (level + 1)); 582 fo = (PetscInt *)malloc(sizeof(PetscInt) * (n + 1)); 583 used = (PetscInt *)malloc(sizeof(PetscInt) * n); 584 585 PCTFS_ivec_zero(dir, level + 1); 586 PCTFS_ivec_zero(nsep, level + 1); 587 PCTFS_ivec_zero(lnsep, level + 1); 588 PCTFS_ivec_set(fo, -1, n + 1); 589 PCTFS_ivec_zero(used, n); 590 591 lhs = (PetscScalar *)malloc(sizeof(PetscScalar) * m); 592 rhs = (PetscScalar *)malloc(sizeof(PetscScalar) * m); 593 594 /* determine the # of unique dof */ 595 PCTFS_rvec_zero(lhs, m); 596 PCTFS_rvec_set(lhs, 1.0, n); 597 PCTFS_gs_gop_hc(PCTFS_gs_handle, lhs, "+\0", level); 598 PetscCall(PetscInfo(0, "done first PCTFS_gs_gop_hc\n")); 599 PCTFS_rvec_zero(rsum, 2); 600 for (i = 0; i < n; i++) { 601 if (lhs[i] != 0.0) { 602 rsum[0] += 1.0 / lhs[i]; 603 rsum[1] += lhs[i]; 604 } 605 if (lhs[i] != 1.0) shared = 1; 606 } 607 608 PCTFS_grop_hc(rsum, rw, 2, op, level); 609 rsum[0] += 0.1; 610 rsum[1] += 0.1; 611 612 xyt_handle->info->n_global = xyt_handle->info->m_global = (PetscInt)rsum[0]; 613 xyt_handle->mvi->n_global = xyt_handle->mvi->m_global = (PetscInt)rsum[0]; 614 615 /* determine separator sets top down */ 616 if (shared) { 617 /* solution is to do as in the symmetric shared case but then */ 618 /* pick the sub-hc with the most free dofs and do a mat-vec */ 619 /* and pick up the responses on the other sub-hc from the */ 620 /* initial separator set obtained from the symm. shared case */ 621 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "shared dof separator determination not ready ... see hmt!!!"); 622 /* [dead code deleted since it is unlikely to be completed] */ 623 } else { 624 for (iptr = fo + n, id = PCTFS_my_id, mask = PCTFS_num_nodes >> 1, edge = level; edge > 0; edge--, mask >>= 1) { 625 /* set rsh of hc, fire, and collect lhs responses */ 626 (id < mask) ? PCTFS_rvec_zero(lhs, m) : PCTFS_rvec_set(lhs, 1.0, m); 627 PCTFS_gs_gop_hc(PCTFS_gs_handle, lhs, "+\0", edge); 628 629 /* set lsh of hc, fire, and collect rhs responses */ 630 (id < mask) ? PCTFS_rvec_set(rhs, 1.0, m) : PCTFS_rvec_zero(rhs, m); 631 PCTFS_gs_gop_hc(PCTFS_gs_handle, rhs, "+\0", edge); 632 633 /* count number of dofs I own that have signal and not in sep set */ 634 for (PCTFS_ivec_zero(sum, 4), ct = i = 0; i < n; i++) { 635 if (!used[i]) { 636 /* number of unmarked dofs on node */ 637 ct++; 638 /* number of dofs to be marked on lhs hc */ 639 if ((id < mask) && (lhs[i] != 0.0)) sum[0]++; 640 /* number of dofs to be marked on rhs hc */ 641 if ((id >= mask) && (rhs[i] != 0.0)) sum[1]++; 642 } 643 } 644 645 /* for the non-symmetric case we need separators of width 2 */ 646 /* so take both sides */ 647 (id < mask) ? (sum[2] = ct) : (sum[3] = ct); 648 PCTFS_giop_hc(sum, w, 4, op, edge); 649 650 ct = 0; 651 if (id < mask) { 652 /* mark dofs I own that have signal and not in sep set */ 653 for (i = 0; i < n; i++) { 654 if ((!used[i]) && (lhs[i] != 0.0)) { 655 ct++; 656 nfo++; 657 *--iptr = local2global[i]; 658 used[i] = edge; 659 } 660 } 661 /* LSH hc summation of ct should be sum[0] */ 662 } else { 663 /* mark dofs I own that have signal and not in sep set */ 664 for (i = 0; i < n; i++) { 665 if ((!used[i]) && (rhs[i] != 0.0)) { 666 ct++; 667 nfo++; 668 *--iptr = local2global[i]; 669 used[i] = edge; 670 } 671 } 672 /* RSH hc summation of ct should be sum[1] */ 673 } 674 675 if (ct > 1) PCTFS_ivec_sort(iptr, ct); 676 lnsep[edge] = ct; 677 nsep[edge] = sum[0] + sum[1]; 678 dir[edge] = BOTH; 679 680 /* LATER or we can recur on these to order seps at this level */ 681 /* do we need full set of separators for this? */ 682 683 /* fold rhs hc into lower */ 684 if (id >= mask) id -= mask; 685 } 686 } 687 688 /* level 0 is on processor case - so mark the remainder */ 689 for (ct = i = 0; i < n; i++) { 690 if (!used[i]) { 691 ct++; 692 nfo++; 693 *--iptr = local2global[i]; 694 used[i] = edge; 695 } 696 } 697 if (ct > 1) PCTFS_ivec_sort(iptr, ct); 698 lnsep[edge] = ct; 699 nsep[edge] = ct; 700 dir[edge] = BOTH; 701 702 xyt_handle->info->nsep = nsep; 703 xyt_handle->info->lnsep = lnsep; 704 xyt_handle->info->fo = fo; 705 xyt_handle->info->nfo = nfo; 706 707 free(dir); 708 free(lhs); 709 free(rhs); 710 free(used); 711 PetscFunctionReturn(0); 712 } 713 714 static mv_info *set_mvi(PetscInt *local2global, PetscInt n, PetscInt m, PetscErrorCode (*matvec)(mv_info *, PetscScalar *, PetscScalar *), void *grid_data) { 715 mv_info *mvi; 716 717 mvi = (mv_info *)malloc(sizeof(mv_info)); 718 mvi->n = n; 719 mvi->m = m; 720 mvi->n_global = -1; 721 mvi->m_global = -1; 722 mvi->local2global = (PetscInt *)malloc((m + 1) * sizeof(PetscInt)); 723 724 PCTFS_ivec_copy(mvi->local2global, local2global, m); 725 mvi->local2global[m] = INT_MAX; 726 mvi->matvec = matvec; 727 mvi->grid_data = grid_data; 728 729 /* set xyt communication handle to perform restricted matvec */ 730 mvi->PCTFS_gs_handle = PCTFS_gs_init(local2global, m, PCTFS_num_nodes); 731 732 return (mvi); 733 } 734 735 static PetscErrorCode do_matvec(mv_info *A, PetscScalar *v, PetscScalar *u) { 736 PetscFunctionBegin; 737 A->matvec((mv_info *)A->grid_data, v, u); 738 PetscFunctionReturn(0); 739 } 740