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