自学教程：C++ vec1函数代码示例

static void emit_points(struct brw_clip_compile *c,			bool do_offset ){   struct brw_compile *p = &c->func;   struct brw_indirect v0 = brw_indirect(0, 0);   struct brw_indirect v0ptr = brw_indirect(2, 0);   brw_MOV(p, c->reg.loopcount, c->reg.nr_verts);   brw_MOV(p, get_addr_reg(v0ptr), brw_address(c->reg.inlist));   brw_DO(p, BRW_EXECUTE_1);   {      brw_MOV(p, get_addr_reg(v0), deref_1uw(v0ptr, 0));      brw_ADD(p, get_addr_reg(v0ptr), get_addr_reg(v0ptr), brw_imm_uw(2));      /* draw if edgeflag != 0        */      brw_CMP(p, 	      vec1(brw_null_reg()), BRW_CONDITIONAL_NZ, 	      deref_1f(v0, brw_varying_to_offset(&c->vue_map,                                                 VARYING_SLOT_EDGE)),	      brw_imm_f(0));      brw_IF(p, BRW_EXECUTE_1);      {	 if (do_offset)	    apply_one_offset(c, v0);	 brw_clip_emit_vue(c, v0, BRW_URB_WRITE_ALLOCATE_COMPLETE,                           (_3DPRIM_POINTLIST << URB_WRITE_PRIM_TYPE_SHIFT)                           | URB_WRITE_PRIM_START | URB_WRITE_PRIM_END);      }      brw_ENDIF(p);      brw_set_conditionalmod(p, BRW_CONDITIONAL_NZ);      brw_ADD(p, c->reg.loopcount, c->reg.loopcount, brw_imm_d(-1));   }   brw_WHILE(p);}

void brw_clip_init_clipmask( struct brw_clip_compile *c ){   struct brw_compile *p = &c->func;   struct brw_reg incoming = get_element_ud(c->reg.R0, 2);      /* Shift so that lowest outcode bit is rightmost:     */   brw_SHR(p, c->reg.planemask, incoming, brw_imm_ud(26));   if (c->key.nr_userclip) {      struct brw_reg tmp = retype(vec1(get_tmp(c)), BRW_REGISTER_TYPE_UD);      /* Rearrange userclip outcodes so that they come directly after       * the fixed plane bits.       */      brw_AND(p, tmp, incoming, brw_imm_ud(0x3f<<14));      brw_SHR(p, tmp, tmp, brw_imm_ud(8));      brw_OR(p, c->reg.planemask, c->reg.planemask, tmp);            release_tmp(c, tmp);   }}

void TestVector() {    SpacePoint vec1(10), vec2, tmp;    //	DataVector<double> vec1(10),vec2,tmp;    vec1 = 5;    //	(SavableClass&)vec2=(SavableClass&)vec1;    vec1[5] = 55;    vec2 = vec1 = 10;    fcout << vec2 << "/n" << vec1 << "/n";    fcout << "vec1 " << vec1 << "/n vec2 " << vec2 << "/n";    tmp = vec1 + vec2;    //	DataVector<double> tmp=vec1+vec2;    fcout << tmp << "/n";    fcout << "sum " << vec1 + vec2 << " /n min " << vec2 - vec1 << " /n mul "          << vec1 * vec2 << "/n";    {        FilterTextOut fo("Test1", DataSource::Memory);        vec2 = 100;        fo << vec2;        FilterTextIn fi("Test1", DataSource::Memory);        fi >> vec1;    }    fcout << "vec1" << vec1 << "/n vec2" << vec2 << "/n";    DataVector<SpacePoint> dat(10);    dat[0].SetDim(10);    dat[0] = 5;    for(int k = 1; k < 9; k++)        dat[k] = dat[0];    fcout << dat << "/n" << (void *)&vec2 << "/n";    void *tmpv = NULL;    FilterTextOut fo("Test1", DataSource::Disk);    fo << (void *)&vec1 << 10;    //fo<<10;    fo.CloseBuf();    FilterTextIn fi("Test1", DataSource::Disk);    int tmp_int;   //fi>>tmp_int;fcout<<tmp_int;    fi >> tmpv >> tmp_int;    fcout << tmpv << tmp_int;}

static void test_scans(unsigned int sz){  viennacl::vector<ScalarType> vec1(sz), vec2(sz);  std::cout << "Initialize vector..." << std::endl;  init_vector(vec1);  // INCLUSIVE SCAN  std::cout << " --- Inclusive scan ---" << std::endl;  std::cout << "Separate vectors: ";  viennacl::linalg::inclusive_scan(vec1, vec2);  test_scan_values(vec1, vec2, true);  std::cout << "In-place: ";  vec2 = vec1;  viennacl::linalg::inclusive_scan(vec2);  test_scan_values(vec1, vec2, true);  std::cout << "Inclusive scan tested successfully!" << std::endl << std::endl;  std::cout << "Initialize vector..." << std::endl;  init_vector(vec1);  // EXCLUSIVE SCAN  std::cout << " --- Exclusive scan ---" << std::endl;  std::cout << "Separate vectors: ";  viennacl::linalg::exclusive_scan(vec1, vec2);  test_scan_values(vec1, vec2, false);  std::cout << "In-place: ";  vec2 = vec1;  viennacl::linalg::exclusive_scan(vec2);  test_scan_values(vec1, vec2, false);  std::cout << "Exclusive scan tested successfully!" << std::endl << std::endl;}

static void emit_points(struct brw_clip_compile *c,			GLboolean do_offset ){   struct brw_compile *p = &c->func;   struct brw_instruction *loop;   struct brw_instruction *draw_point;   struct brw_indirect v0 = brw_indirect(0, 0);   struct brw_indirect v0ptr = brw_indirect(2, 0);   brw_MOV(p, c->reg.loopcount, c->reg.nr_verts);   brw_MOV(p, get_addr_reg(v0ptr), brw_address(c->reg.inlist));   loop = brw_DO(p, BRW_EXECUTE_1);   {      brw_MOV(p, get_addr_reg(v0), deref_1uw(v0ptr, 0));      brw_ADD(p, get_addr_reg(v0ptr), get_addr_reg(v0ptr), brw_imm_uw(2));      /* draw if edgeflag != 0        */      brw_CMP(p, 	      vec1(brw_null_reg()), BRW_CONDITIONAL_NZ, 	      deref_1f(v0, c->offset[VERT_RESULT_EDGE]),	      brw_imm_f(0));      draw_point = brw_IF(p, BRW_EXECUTE_1);      {	 if (do_offset)	    apply_one_offset(c, v0);	 brw_clip_emit_vue(c, v0, 1, 0, (_3DPRIM_POINTLIST << 2) | R02_PRIM_START | R02_PRIM_END);      }      brw_ENDIF(p, draw_point);      brw_set_conditionalmod(p, BRW_CONDITIONAL_NZ);      brw_ADD(p, c->reg.loopcount, c->reg.loopcount, brw_imm_d(-1));   }   brw_WHILE(p, loop);}

void brw_clip_tri_init_vertices( struct brw_clip_compile *c ){   struct brw_compile *p = &c->func;   struct brw_reg tmp0 = c->reg.loopcount; /* handy temporary */   struct brw_instruction *is_rev;   /* Initial list of indices for incoming vertexes:    */   brw_AND(p, tmp0, get_element_ud(c->reg.R0, 2), brw_imm_ud(PRIM_MASK));    brw_CMP(p, 	   vec1(brw_null_reg()), 	   BRW_CONDITIONAL_EQ, 	   tmp0,	   brw_imm_ud(_3DPRIM_TRISTRIP_REVERSE));   /* XXX: Is there an easier way to do this?  Need to reverse every    * second tristrip element:  Can ignore sometimes?    */   is_rev = brw_IF(p, BRW_EXECUTE_1);   {         brw_MOV(p, get_element(c->reg.inlist, 0),  brw_address(c->reg.vertex[1]) );      brw_MOV(p, get_element(c->reg.inlist, 1),  brw_address(c->reg.vertex[0]) );      if (c->need_direction)	 brw_MOV(p, c->reg.dir, brw_imm_f(-1));   }   is_rev = brw_ELSE(p, is_rev);   {      brw_MOV(p, get_element(c->reg.inlist, 0),  brw_address(c->reg.vertex[0]) );      brw_MOV(p, get_element(c->reg.inlist, 1),  brw_address(c->reg.vertex[1]) );      if (c->need_direction)	 brw_MOV(p, c->reg.dir, brw_imm_f(1));   }   brw_ENDIF(p, is_rev);   brw_MOV(p, get_element(c->reg.inlist, 2),  brw_address(c->reg.vertex[2]) );   brw_MOV(p, brw_vec8_grf(c->reg.outlist.nr, 0), brw_imm_f(0));   brw_MOV(p, c->reg.nr_verts, brw_imm_ud(3));}

	void namevalue_object_t::test<20>()	{		LLNameValue nValue1(" SecondLife S32 RW SIM 22222");		LLNameValue nValue2(" Virtual S32 RW SIM 33333");		LLNameValue nValue3(" SecondLife S32"); 		nValue3 = nValue1 % nValue2;		ensure_equals("1:operator% failed",*nValue3.getS32(),22222);		LLNameValue nValue4(" SecondLife U32 RW SIM 3");		LLNameValue nValue5(" SecondLife S32 RW SIM 2");		LLNameValue nValue6(" SecondLife S32");		nValue6 = nValue4 % nValue5;		ensure_equals("2:operator% failed",*nValue6.getS32(),1);		LLNameValue nValue10(" SecondLife VEC3 RW SIM <4, 5, 6>");		LLNameValue nValue11(" SecondLife VEC3 RW SIM <1, 2, 3>");		LLNameValue nValue12(" SecondLife VEC3");		LLVector3 vec1(4,5,6);		LLVector3 vec2(1,2,3);		LLVector3 vec3(vec1 % vec2);		nValue12 = nValue10 % nValue11;		ensure_equals("5:operator% failed",*nValue12.getVec3(), vec3);	}

// see http://llvm.org/docs/LibFuzzer.htmlextern "C" int LLVMFuzzerTestOneInput(const uint8_t* data, size_t size){    try    {        // step 1: parse input        std::vector<uint8_t> vec1(data, data + size);        json j1 = json::from_cbor(vec1);        try        {            // step 2: round trip            std::vector<uint8_t> vec2 = json::to_cbor(j1);            // parse serialization            json j2 = json::from_cbor(vec2);            // serializations must match            assert(json::to_cbor(j2) == vec2);        }        catch (const json::parse_error&)        {            // parsing a CBOR serialization must not fail            assert(false);        }    }    catch (const json::parse_error&)    {        // parse errors are ok, because input may be random bytes    }    catch (const json::type_error&)    {        // type errors can occur during parsing, too    }    // return 0 - non-zero return values are reserved for future use    return 0;}

TEST_F(UUIDTest, test6){    std::list<UUID> vec1(20);    const fs::path tmpfile(FileUtils::create_temp_file(FileUtils::temp_path(),                                                       "serialization.txt"));    ALWAYS_CLEANUP_FILE(tmpfile);    {        fs::ofstream ofs(tmpfile);        boost::archive::xml_oarchive oa(ofs);        oa << BOOST_SERIALIZATION_NVP(vec1);    }    std::list<UUID> vec2;    fs::ifstream ifs(tmpfile);    boost::archive::xml_iarchive ia(ifs);    ia >> BOOST_SERIALIZATION_NVP(vec2);    ASSERT_TRUE(vec1 == vec2);}

void TestVectorManipulation(size_t n){    typedef typename Vector::value_type T;    thrust::host_vector<T> src = unittest::random_samples<T>(n);    ASSERT_EQUAL(src.size(), n);    // basic initialization    Vector test0(n);    Vector test1(n, (T) 3);    ASSERT_EQUAL(test0.size(), n);    ASSERT_EQUAL(test1.size(), n);    ASSERT_EQUAL((test1 == std::vector<T>(n, (T) 3)), true);       // initializing from other vector    std::vector<T> stl_vector(src.begin(), src.end());    Vector cpy0 = src;    Vector cpy1(stl_vector);    Vector cpy2(stl_vector.begin(), stl_vector.end());    ASSERT_EQUAL(cpy0, src);    ASSERT_EQUAL(cpy1, src);    ASSERT_EQUAL(cpy2, src);    // resizing    Vector vec1(src);    vec1.resize(n + 3);    ASSERT_EQUAL(vec1.size(), n + 3);    vec1.resize(n);    ASSERT_EQUAL(vec1.size(), n);    ASSERT_EQUAL(vec1, src);         vec1.resize(n + 20, (T) 11);    Vector tail(vec1.begin() + n, vec1.end());    ASSERT_EQUAL( (tail == std::vector<T>(20, (T) 11)), true);    vec1.resize(0);    ASSERT_EQUAL(vec1.size(), 0);    ASSERT_EQUAL(vec1.empty(), true);    vec1.resize(10);    ASSERT_EQUAL(vec1.size(), 10);    vec1.clear();    ASSERT_EQUAL(vec1.size(), 0);    vec1.resize(5);    ASSERT_EQUAL(vec1.size(), 5);    // push_back    Vector vec2;    for(size_t i = 0; i < 10; ++i){        ASSERT_EQUAL(vec2.size(), i);        vec2.push_back( (T) i );        ASSERT_EQUAL(vec2.size(), i + 1);        for(size_t j = 0; j <= i; j++)            ASSERT_EQUAL(vec2[j],     j);        ASSERT_EQUAL(vec2.back(), i);    }    // pop_back    for(size_t i = 10; i > 0; --i){        ASSERT_EQUAL(vec2.size(), i);        ASSERT_EQUAL(vec2.back(), i-1);        vec2.pop_back();        ASSERT_EQUAL(vec2.size(), i-1);        for(size_t j = 0; j < i; j++)            ASSERT_EQUAL(vec2[j], j);    }    //TODO test swap, erase(pos), erase(begin, end)}

int main(){  typedef float       ScalarType;  //  // Initialize OpenCL vectors:  //  unsigned int vector_size = 10;  viennacl::scalar<ScalarType>  s = 1.0; //dummy  viennacl::vector<ScalarType>  vec1(vector_size);  viennacl::vector<ScalarType>  vec2(vector_size);  viennacl::vector<ScalarType>  result_mul(vector_size);  viennacl::vector<ScalarType>  result_div(vector_size);  //  // fill the operands vec1 and vec2:  //  for (unsigned int i=0; i<vector_size; ++i)  {    vec1[i] = static_cast<ScalarType>(i);    vec2[i] = static_cast<ScalarType>(vector_size-i);  }  //  // Set up the OpenCL program given in my_compute_kernel:  // A program is one compilation unit and can hold many different compute kernels.  //  viennacl::ocl::program & my_prog = viennacl::ocl::current_context().add_program(my_compute_program, "my_compute_program");  my_prog.add_kernel("elementwise_prod");  //register elementwise product kernel  my_prog.add_kernel("elementwise_div");   //register elementwise division kernel    //  // Now we can get the kernels from the program 'my_program'.  // (Note that first all kernels need to be registered via add_kernel() before get_kernel() can be called,  // otherwise existing references might be invalidated)  //  viennacl::ocl::kernel & my_kernel_mul = my_prog.get_kernel("elementwise_prod");  viennacl::ocl::kernel & my_kernel_div = my_prog.get_kernel("elementwise_div");    //  // Launch the kernel with 'vector_size' threads in one work group  // Note that size_t might differ between host and device. Thus, a cast to cl_uint is necessary for the forth argument.  //  viennacl::ocl::enqueue(my_kernel_mul(vec1, vec2, result_mul, static_cast<cl_uint>(vec1.size())));    viennacl::ocl::enqueue(my_kernel_div(vec1, vec2, result_div, static_cast<cl_uint>(vec1.size())));    //  // Print the result:  //  std::cout << "        vec1: " << vec1 << std::endl;  std::cout << "        vec2: " << vec2 << std::endl;  std::cout << "vec1 .* vec2: " << result_mul << std::endl;  std::cout << "vec1 /* vec2: " << result_div << std::endl;  std::cout << "norm_2(vec1 .* vec2): " << viennacl::linalg::norm_2(result_mul) << std::endl;  std::cout << "norm_2(vec1 /* vec2): " << viennacl::linalg::norm_2(result_div) << std::endl;    //  //  That's it.  //  std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;    return 0;}

//.........这里部分代码省略.........	 if (c->dispatch_width == 16) {	    brw_set_compression_control(p, BRW_COMPRESSION_2NDHALF);	    brw_MOV(p,		    brw_message_reg(nr + channel + 4),		    sechalf(arg0[channel]));	 }      }   }   brw_set_saturate(p, 0);   /* skip over the regs populated above:    */   if (c->dispatch_width == 16)      nr += 8;   else      nr += 4;   brw_pop_insn_state(p);   if (c->source_depth_to_render_target)   {      if (c->computes_depth)	 brw_MOV(p, brw_message_reg(nr), arg2[2]);      else 	 brw_MOV(p, brw_message_reg(nr), arg1[1]); /* ? */      nr += 2;   }   if (c->dest_depth_reg)   {      GLuint comp = c->dest_depth_reg / 2;      GLuint off = c->dest_depth_reg % 2;      if (off != 0) {         brw_push_insn_state(p);         brw_set_compression_control(p, BRW_COMPRESSION_NONE);         brw_MOV(p, brw_message_reg(nr), offset(arg1[comp],1));         /* 2nd half? */         brw_MOV(p, brw_message_reg(nr+1), arg1[comp+1]);         brw_pop_insn_state(p);      }      else {         brw_MOV(p, brw_message_reg(nr), arg1[comp]);      }      nr += 2;   }   if (intel->gen >= 6) {      /* Load the message header.  There's no implied move from src0       * to the base mrf on gen6.       */      brw_push_insn_state(p);      brw_set_mask_control(p, BRW_MASK_DISABLE);      brw_MOV(p, retype(brw_message_reg(0), BRW_REGISTER_TYPE_UD),	      retype(brw_vec8_grf(0, 0), BRW_REGISTER_TYPE_UD));      brw_pop_insn_state(p);      if (target != 0) {	 brw_MOV(p, retype(brw_vec1_reg(BRW_MESSAGE_REGISTER_FILE,					0,					2), BRW_REGISTER_TYPE_UD),		 brw_imm_ud(target));      }   }   if (!c->runtime_check_aads_emit) {      if (c->aa_dest_stencil_reg)	 emit_aa(c, arg1, 2);      fire_fb_write(c, 0, nr, target, eot);   }   else {      struct brw_reg v1_null_ud = vec1(retype(brw_null_reg(), BRW_REGISTER_TYPE_UD));      struct brw_reg ip = brw_ip_reg();      struct brw_instruction *jmp;            brw_set_compression_control(p, BRW_COMPRESSION_NONE);      brw_set_conditionalmod(p, BRW_CONDITIONAL_Z);      brw_AND(p, 	      v1_null_ud, 	      get_element_ud(brw_vec8_grf(1,0), 6), 	      brw_imm_ud(1<<26));       jmp = brw_JMPI(p, ip, ip, brw_imm_w(0));      {	 emit_aa(c, arg1, 2);	 fire_fb_write(c, 0, nr, target, eot);	 /* note - thread killed in subroutine */      }      brw_land_fwd_jump(p, jmp);      /* ELSE: Shuffle up one register to fill in the hole left for AA:       */      fire_fb_write(c, 1, nr-1, target, eot);   }}

void nuiGLDrawContext::DrawGradient(const nuiGradient& rGradient, const nuiRect& rEnclosingRect, nuiSize x1, nuiSize y1, nuiSize x2, nuiSize y2){  nglVector2f vec(x2 - x1, y2 - y1);  nglVector2f para(-vec[1], vec[0]);  nglVector2f vec1(vec);  nglVector2f para1(para);  vec1.Normalize();  para1.Normalize();  // What Quadrant are we in?:  //         |  //     a   |   b  //         |  //  ----------------  //         |  //     c   |   d  //         |  float xa, xb, xc, xd;  float ya, yb, yc, yd;  float x, y;  float xp, yp;  float xx, yy;  float xxp, yyp;  xa = xc = rEnclosingRect.Left();  xb = xd = rEnclosingRect.Right();  ya = yb = rEnclosingRect.Top();  yc = yd = rEnclosingRect.Bottom();  if (x1 < x2)  {    // Go from a to d or c to b    if (y1 == y2)    {      x  = xa; y  = ya;      xp = xc; yp = yc;      xx = xd; yy = yd;      xxp= xb; yyp= yb;    }    else if (y1 < y2)    {      // a to d      IntersectLines(xa,ya, para1[0], para1[1], xb, yb, vec1[0], vec1[1], x, y);      IntersectLines(xa,ya, para1[0], para1[1], xc, yc, vec1[0], vec1[1], xp, yp);      IntersectLines(xd,yd, para1[0], para1[1], xc, yc, vec1[0], vec1[1], xx, yy);      IntersectLines(xd,yd, para1[0], para1[1], xb, yb, vec1[0], vec1[1], xxp, yyp);    }    else    {      // c to d      IntersectLines(xc,yc, para1[0], para1[1], xa, ya, vec1[0], vec1[1], x, y);      IntersectLines(xc,yc, para1[0], para1[1], xd, yd, vec1[0], vec1[1], xp, yp);      IntersectLines(xb,yb, para1[0], para1[1], xd, yd, vec1[0], vec1[1], xx, yy);      IntersectLines(xb,yb, para1[0], para1[1], xa, ya, vec1[0], vec1[1], xxp, yyp);    }  }  else  {    if (y1 == y2)    {      x  = xd; y  = yd;      xp = xb; yp = yb;      xx = xa; yy = ya;      xxp= xc; yyp= yc;    }    else if (y1 < y2)    {      // b to c      IntersectLines(xb,yb, para1[0], para1[1], xd, yd, vec1[0], vec1[1], x, y);      IntersectLines(xb,yb, para1[0], para1[1], xa, ya, vec1[0], vec1[1], xp, yp);      IntersectLines(xc,yc, para1[0], para1[1], xa, ya, vec1[0], vec1[1], xx, yy);      IntersectLines(xc,yc, para1[0], para1[1], xd, yd, vec1[0], vec1[1], xxp, yyp);    }    else    {      // d to a      IntersectLines(xd,yd, para1[0], para1[1], xc, yc, vec1[0], vec1[1], x, y);      IntersectLines(xd,yd, para1[0], para1[1], xb, yb, vec1[0], vec1[1], xp, yp);      IntersectLines(xa,ya, para1[0], para1[1], xb, yb, vec1[0], vec1[1], xx, yy);      IntersectLines(xa,ya, para1[0], para1[1], xc, yc, vec1[0], vec1[1], xxp, yyp);    }  }  float startx,starty;  float startxp,startyp;  float stopx,stopy;  float stopxp,stopyp;  if (y1 != y2)  {    IntersectLines(x1, y1, para1[0], para1[1], x,  y,  vec1[0], vec1[1], startx,  starty);    IntersectLines(x1, y1, para1[0], para1[1], xp, yp, vec1[0], vec1[1], startxp, startyp);    IntersectLines(x2, y2, para1[0], para1[1], x,  y,  vec1[0], vec1[1], stopx,   stopy);    IntersectLines(x2, y2, para1[0], para1[1], xp, yp, vec1[0], vec1[1], stopxp,  stopyp);  }  else  {    startx  = x1; starty  = y;    startxp = x1; startyp = yp;    stopx   = x2; stopy   = y;//.........这里部分代码省略.........

/* Interpolate between two vertices and put the result into a0.0. * Increment a0.0 accordingly. * * Beware that dest_ptr can be equal to v0_ptr! */void brw_clip_interp_vertex( struct brw_clip_compile *c,			     struct brw_indirect dest_ptr,			     struct brw_indirect v0_ptr, /* from */			     struct brw_indirect v1_ptr, /* to */			     struct brw_reg t0,			     bool force_edgeflag){   struct brw_codegen *p = &c->func;   struct brw_reg t_nopersp, v0_ndc_copy;   GLuint slot;   /* Just copy the vertex header:    */   /*    * After CLIP stage, only first 256 bits of the VUE are read    * back on Ironlake, so needn't change it    */   brw_copy_indirect_to_indirect(p, dest_ptr, v0_ptr, 1);   /* First handle the 3D and NDC interpolation, in case we    * need noperspective interpolation. Doing it early has no    * performance impact in any case.    */   /* Take a copy of the v0 NDC coordinates, in case dest == v0. */   if (c->has_noperspective_shading) {      GLuint offset = brw_varying_to_offset(&c->vue_map,                                                 BRW_VARYING_SLOT_NDC);      v0_ndc_copy = get_tmp(c);      brw_MOV(p, v0_ndc_copy, deref_4f(v0_ptr, offset));   }   /* Compute the new 3D position    *    * dest_hpos = v0_hpos * (1 - t0) + v1_hpos * t0    */   {      GLuint delta = brw_varying_to_offset(&c->vue_map, VARYING_SLOT_POS);      struct brw_reg tmp = get_tmp(c);      brw_MUL(p, vec4(brw_null_reg()), deref_4f(v1_ptr, delta), t0);      brw_MAC(p, tmp, negate(deref_4f(v0_ptr, delta)), t0);      brw_ADD(p, deref_4f(dest_ptr, delta), deref_4f(v0_ptr, delta), tmp);      release_tmp(c, tmp);   }   /* Recreate the projected (NDC) coordinate in the new vertex header */   brw_clip_project_vertex(c, dest_ptr);   /* If we have noperspective attributes,    * we need to compute the screen-space t    */   if (c->has_noperspective_shading) {      GLuint delta = brw_varying_to_offset(&c->vue_map,                                                BRW_VARYING_SLOT_NDC);      struct brw_reg tmp = get_tmp(c);      t_nopersp = get_tmp(c);      /* t_nopersp = vec4(v1.xy, dest.xy) */      brw_MOV(p, t_nopersp, deref_4f(v1_ptr, delta));      brw_MOV(p, tmp, deref_4f(dest_ptr, delta));      brw_set_default_access_mode(p, BRW_ALIGN_16);      brw_MOV(p,              brw_writemask(t_nopersp, WRITEMASK_ZW),              brw_swizzle(tmp, 0, 1, 0, 1));      /* t_nopersp = vec4(v1.xy, dest.xy) - v0.xyxy */      brw_ADD(p, t_nopersp, t_nopersp,              negate(brw_swizzle(v0_ndc_copy, 0, 1, 0, 1)));      /* Add the absolute values of the X and Y deltas so that if       * the points aren't in the same place on the screen we get       * nonzero values to divide.       *       * After that, we have vert1 - vert0 in t_nopersp.x and       * vertnew - vert0 in t_nopersp.y       *       * t_nopersp = vec2(|v1.x  -v0.x| + |v1.y  -v0.y|,       *                  |dest.x-v0.x| + |dest.y-v0.y|)       */      brw_ADD(p,              brw_writemask(t_nopersp, WRITEMASK_XY),              brw_abs(brw_swizzle(t_nopersp, 0, 2, 0, 0)),              brw_abs(brw_swizzle(t_nopersp, 1, 3, 0, 0)));      brw_set_default_access_mode(p, BRW_ALIGN_1);      /* If the points are in the same place, just substitute a       * value to avoid divide-by-zero       */      brw_CMP(p, vec1(brw_null_reg()), BRW_CONDITIONAL_EQ,              vec1(t_nopersp),              brw_imm_f(0));      brw_IF(p, BRW_EXECUTE_1);      brw_MOV(p, t_nopersp, brw_imm_vf4(brw_float_to_vf(1.0),                                        brw_float_to_vf(0.0),//.........这里部分代码省略.........

/*** Example using STL adaptors*/void STLMoblet::STL_adaptors(){	LOG("/n");	LOG("========================= STL adaptors ==========================================================================");	LOG("/**");	LOG("* Function object adaptors are used to create a function object from another function object."				);	LOG("* The created function object, is not the same as the original functor, but is adapted to a certain need."	);	LOG("* For example if we have a function object like std::less, and we want to compare all the elements"		);	LOG("* in a range against 10, then be can use an STL adaptor (bind2nd), that bounds the second argument"		);	LOG("* of std::less to 10. Then we can use std::less with algorithms like remove_copy_if,"						);	LOG("* that need a unary predicate."																			);	LOG("*"																											);	LOG("* STL provides two functor adaptors: bind1st and bind2nd."													);	LOG("*"																											);	LOG("* 	bind1st: constructs an unary function object from a binary function object, by binding"					);	LOG("* 	the first argument to a fixed value."																	);	LOG("*"																											);	LOG("* 	bind1st template function is defined like this:"														);	LOG("*"																											);	LOG("*   binder1st<SomeFunctor> bind1st (const SomeFunctor& fun, const T& fixedValue"							);	LOG("*	{"																										);	LOG("*		return binder1st<SomeFunctor>(fun, x);"																);	LOG("*	}"																										);	LOG("*"																											);	LOG("*	bind1st returns an binder1st object, which is actually a functor that forwards the function calls"		);	LOG("*	to the /"fun/" argument it takes as a parameter, when constructed:"										);	LOG("*"																											);	LOG("*	template <class Functor> class binder1st {"																);	LOG("*"																											);	LOG("*		binder1st(const Functor &fun, Functor::first_argument_type &fixed)"									);	LOG("*		{"																									);	LOG("*			mFun = fun;"																					);	LOG("*			mFixedValue = fixed;"																			);	LOG("*		}"																									);	LOG("*"																											);	LOG("*		Functor::result_type operator()(Functor::second_argument_type &someValue){"							);	LOG("*			return mFun(mFixedValue, someValue);"															);	LOG("*		}"																									);	LOG("*	};"																										);	LOG("*"																											);	LOG("*	bind2nd is implemented in a similar way, but instead of binding the first argument,"					);	LOG("*	bind2nd it will bind the second one to a fixed value."													);	LOG("*"																											);	LOG("*	bind1st and bind2nd function templates are defined in the <functional> header."							);	LOG("*/");	LOG("/n");	LOG("				Example using adaptors 							  ");	LOG("/n																				   ");	LOG("/**"																				);	LOG("*  bind1st function template: constructs an unary function object from a"			);	LOG("*  binary function object, by binding the first parameter to a certain value."		);	LOG("/n    */"																			);	log_to_console("/n     Example using std::bind1st:/n");	int array[] = { 1, -99, 2, -100 };	int arraySize = sizeof(array)/sizeof(array[0]);	std::vector<int> vec1(array, array + arraySize);	log_to_console(vec1, "vec1 contains: ");	LOG("/n"																		);	LOG("/**std::remove_if calls std::less(0, element). If less(0,element) returns"	);	LOG("*  true => removes that element."											);	LOG("*  less(0,element) is equivalent to 0<element."							);	LOG("*/"																		);	LOG("/n"																		);	TRACE(std::vector<int>::iterator newEnd = std::remove_if(vec1.begin(),			vec1.end(), bind1st(std::less<int>(), 1)));	log_to_console("vec1 after calling std::remove_if(vec1.begin(), vec1.end(), "			"bind1st(std::less<int>(), 0)): ");	for(std::vector<int>::iterator it = vec1.begin(); it != newEnd; ++it)	{		log_to_console(*it);	}	LOG("/n"															      	  );	LOG("/**"																  	  );	LOG("*  bind2nd function template: constructs an unary function object from a");	LOG("*  binary function object, by binding the second parameter to a certain" );	LOG("*  value."															 	  );	LOG("*/"																	  );	log_to_console("/n     Example using std::bind2nd:/n");	std::vector<int> vec2(array, array + arraySize);	log_to_console(vec2, "vec2 contains: ");	LOG("/n"																	  );	LOG("/**  std::remove_if calls std::greater(element, 0)."					  );//.........这里部分代码省略.........

void planning_environment::setMarkerShapeFromShape(const shapes::Shape *obj, visualization_msgs::Marker &mk, double padding){  switch (obj->type)  {  case shapes::SPHERE:    mk.type = visualization_msgs::Marker::SPHERE;    mk.scale.x = mk.scale.y = mk.scale.z = static_cast<const shapes::Sphere*>(obj)->radius * 2.0 + padding;    break;      case shapes::BOX:    mk.type = visualization_msgs::Marker::CUBE;    {      const double *size = static_cast<const shapes::Box*>(obj)->size;      mk.scale.x = size[0] + padding*2.0;      mk.scale.y = size[1] + padding*2.0;      mk.scale.z = size[2] + padding*2.0;    }    break;      case shapes::CYLINDER:    mk.type = visualization_msgs::Marker::CYLINDER;    mk.scale.x = static_cast<const shapes::Cylinder*>(obj)->radius * 2.0 + padding;    mk.scale.y = mk.scale.x;    mk.scale.z = static_cast<const shapes::Cylinder*>(obj)->length + padding*2.0;    break;      case shapes::MESH:    mk.type = visualization_msgs::Marker::LINE_LIST;    mk.scale.x = mk.scale.y = mk.scale.z = 0.001;    {	         const shapes::Mesh *mesh = static_cast<const shapes::Mesh*>(obj);      double* vertices = new double[mesh->vertexCount * 3];      double sx = 0.0, sy = 0.0, sz = 0.0;      for(unsigned int i = 0; i < mesh->vertexCount; ++i) {        unsigned int i3 = i * 3;        vertices[i3] = mesh->vertices[i3];        vertices[i3 + 1] = mesh->vertices[i3 + 1];        vertices[i3 + 2] = mesh->vertices[i3 + 2];        sx += vertices[i3];        sy += vertices[i3 + 1];        sz += vertices[i3 + 2];      }      // the center of the mesh      sx /= (double)mesh->vertexCount;      sy /= (double)mesh->vertexCount;      sz /= (double)mesh->vertexCount;            for (unsigned int i = 0 ; i < mesh->vertexCount ; ++i)      {        unsigned int i3 = i * 3;	        // vector from center to the vertex        double dx = vertices[i3] - sx;        double dy = vertices[i3 + 1] - sy;        double dz = vertices[i3 + 2] - sz;	        // length of vector        //double norm = sqrt(dx * dx + dy * dy + dz * dz);	        double ndx = ((dx > 0) ? dx+padding : dx-padding);        double ndy = ((dy > 0) ? dy+padding : dy-padding);        double ndz = ((dz > 0) ? dz+padding : dz-padding);                // the new distance of the vertex from the center        //double fact = scale + padding/norm;        vertices[i3] = sx + ndx; //dx * fact;        vertices[i3 + 1] = sy + ndy; //dy * fact;        vertices[i3 + 2] = sz + ndz; //dz * fact;		          }            tf::Transform trans;      tf::poseMsgToTF(mk.pose, trans);      for (unsigned int j = 0 ; j < mesh->triangleCount; ++j) {        unsigned int t1ind = mesh->triangles[3*j];        unsigned int t2ind = mesh->triangles[3*j + 1];        unsigned int t3ind = mesh->triangles[3*j + 2];                tf::Vector3 vec1(vertices[t1ind*3],                       vertices[t1ind*3+1],                       vertices[t1ind*3+2]);                tf::Vector3 vec2(vertices[t2ind*3],                       vertices[t2ind*3+1],                       vertices[t2ind*3+2]);                tf::Vector3 vec3(vertices[t3ind*3],                       vertices[t3ind*3+1],                       vertices[t3ind*3+2]);                //vec1 = trans*vec1;        //vec2 = trans*vec2;        //vec3 = trans*vec3;                geometry_msgs::Point pt1;        pt1.x = vec1.x();        pt1.y = vec1.y();        pt1.z = vec1.z();        geometry_msgs::Point pt2;//.........这里部分代码省略.........

/** * Generate the geometry shader program used on Gen6 to perform stream output * (transform feedback). */voidgen6_sol_program(struct brw_gs_compile *c, struct brw_gs_prog_key *key,	         unsigned num_verts, bool check_edge_flags){   struct brw_compile *p = &c->func;   c->prog_data.svbi_postincrement_value = num_verts;   brw_gs_alloc_regs(c, num_verts, true);   brw_gs_initialize_header(c);   if (key->num_transform_feedback_bindings > 0) {      unsigned vertex, binding;      struct brw_reg destination_indices_uw =         vec8(retype(c->reg.destination_indices, BRW_REGISTER_TYPE_UW));      /* Note: since we use the binding table to keep track of buffer offsets       * and stride, the GS doesn't need to keep track of a separate pointer       * into each buffer; it uses a single pointer which increments by 1 for       * each vertex.  So we use SVBI0 for this pointer, regardless of whether       * transform feedback is in interleaved or separate attribs mode.       *       * Make sure that the buffers have enough room for all the vertices.       */      brw_ADD(p, get_element_ud(c->reg.temp, 0),	         get_element_ud(c->reg.SVBI, 0), brw_imm_ud(num_verts));      brw_CMP(p, vec1(brw_null_reg()), BRW_CONDITIONAL_LE,	         get_element_ud(c->reg.temp, 0),	         get_element_ud(c->reg.SVBI, 4));      brw_IF(p, BRW_EXECUTE_1);      /* Compute the destination indices to write to.  Usually we use SVBI[0]       * + (0, 1, 2).  However, for odd-numbered triangles in tristrips, the       * vertices come down the pipeline in reversed winding order, so we need       * to flip the order when writing to the transform feedback buffer.  To       * ensure that flatshading accuracy is preserved, we need to write them       * in order SVBI[0] + (0, 2, 1) if we're using the first provoking       * vertex convention, and in order SVBI[0] + (1, 0, 2) if we're using       * the last provoking vertex convention.       *       * Note: since brw_imm_v can only be used in instructions in       * packed-word execution mode, and SVBI is a double-word, we need to       * first move the appropriate immediate constant ((0, 1, 2), (0, 2, 1),       * or (1, 0, 2)) to the destination_indices register, and then add SVBI       * using a separate instruction.  Also, since the immediate constant is       * expressed as packed words, and we need to load double-words into       * destination_indices, we need to intersperse zeros to fill the upper       * halves of each double-word.       */      brw_MOV(p, destination_indices_uw,              brw_imm_v(0x00020100)); /* (0, 1, 2) */      if (num_verts == 3) {         /* Get primitive type into temp register. */         brw_AND(p, get_element_ud(c->reg.temp, 0),                 get_element_ud(c->reg.R0, 2), brw_imm_ud(0x1f));         /* Test if primitive type is TRISTRIP_REVERSE.  We need to do this as          * an 8-wide comparison so that the conditional MOV that follows          * moves all 8 words correctly.          */         brw_CMP(p, vec8(brw_null_reg()), BRW_CONDITIONAL_EQ,                 get_element_ud(c->reg.temp, 0),                 brw_imm_ud(_3DPRIM_TRISTRIP_REVERSE));         /* If so, then overwrite destination_indices_uw with the appropriate          * reordering.          */         brw_MOV(p, destination_indices_uw,                 brw_imm_v(key->pv_first ? 0x00010200    /* (0, 2, 1) */                                         : 0x00020001)); /* (1, 0, 2) */         brw_set_predicate_control(p, BRW_PREDICATE_NONE);      }      brw_ADD(p, c->reg.destination_indices,              c->reg.destination_indices, get_element_ud(c->reg.SVBI, 0));      /* For each vertex, generate code to output each varying using the       * appropriate binding table entry.       */      for (vertex = 0; vertex < num_verts; ++vertex) {         /* Set up the correct destination index for this vertex */         brw_MOV(p, get_element_ud(c->reg.header, 5),                 get_element_ud(c->reg.destination_indices, vertex));         for (binding = 0; binding < key->num_transform_feedback_bindings;              ++binding) {            unsigned char varying =               key->transform_feedback_bindings[binding];            unsigned char slot = c->vue_map.varying_to_slot[varying];            /* From the Sandybridge PRM, Volume 2, Part 1, Section 4.5.1:             *             *   "Prior to End of Thread with a URB_WRITE, the kernel must             *   ensure that all writes are complete by sending the final             *   write as a committed write."             */            bool final_write =               binding == key->num_transform_feedback_bindings - 1 &&               vertex == num_verts - 1;//.........这里部分代码省略.........

/* Line clipping, more or less following the following algorithm: * *  for (p=0;p<MAX_PLANES;p++) { *     if (clipmask & (1 << p)) { *        GLfloat dp0 = DOTPROD( vtx0, plane[p] ); *        GLfloat dp1 = DOTPROD( vtx1, plane[p] ); * *        if (dp1 < 0.0f) { *           GLfloat t = dp1 / (dp1 - dp0); *           if (t > t1) t1 = t; *        } else { *           GLfloat t = dp0 / (dp0 - dp1); *           if (t > t0) t0 = t; *        } * *        if (t0 + t1 >= 1.0) *           return; *     } *  } * *  interp( ctx, newvtx0, vtx0, vtx1, t0 ); *  interp( ctx, newvtx1, vtx1, vtx0, t1 ); * */static void clip_and_emit_line( struct brw_clip_compile *c ){   struct brw_codegen *p = &c->func;   struct brw_indirect vtx0     = brw_indirect(0, 0);   struct brw_indirect vtx1      = brw_indirect(1, 0);   struct brw_indirect newvtx0   = brw_indirect(2, 0);   struct brw_indirect newvtx1   = brw_indirect(3, 0);   struct brw_indirect plane_ptr = brw_indirect(4, 0);   struct brw_reg v1_null_ud = retype(vec1(brw_null_reg()), BRW_REGISTER_TYPE_UD);   GLuint hpos_offset = brw_varying_to_offset(&c->vue_map, VARYING_SLOT_POS);   GLint clipdist0_offset = c->key.nr_userclip      ? brw_varying_to_offset(&c->vue_map, VARYING_SLOT_CLIP_DIST0)      : 0;   brw_MOV(p, get_addr_reg(vtx0),      brw_address(c->reg.vertex[0]));   brw_MOV(p, get_addr_reg(vtx1),      brw_address(c->reg.vertex[1]));   brw_MOV(p, get_addr_reg(newvtx0),   brw_address(c->reg.vertex[2]));   brw_MOV(p, get_addr_reg(newvtx1),   brw_address(c->reg.vertex[3]));   brw_MOV(p, get_addr_reg(plane_ptr), brw_clip_plane0_address(c));   /* Note: init t0, t1 together:    */   brw_MOV(p, vec2(c->reg.t0), brw_imm_f(0));   brw_clip_init_planes(c);   brw_clip_init_clipmask(c);   /* -ve rhw workaround */   if (p->devinfo->has_negative_rhw_bug) {      brw_AND(p, brw_null_reg(), get_element_ud(c->reg.R0, 2),              brw_imm_ud(1<<20));      brw_inst_set_cond_modifier(p->devinfo, brw_last_inst, BRW_CONDITIONAL_NZ);      brw_OR(p, c->reg.planemask, c->reg.planemask, brw_imm_ud(0x3f));      brw_inst_set_pred_control(p->devinfo, brw_last_inst, BRW_PREDICATE_NORMAL);   }   /* Set the initial vertex source mask: The first 6 planes are the bounds    * of the view volume; the next 8 planes are the user clipping planes.    */   brw_MOV(p, c->reg.vertex_src_mask, brw_imm_ud(0x3fc0));   /* Set the initial clipdistance offset to be 6 floats before gl_ClipDistance[0].    * We'll increment 6 times before we start hitting actual user clipping. */   brw_MOV(p, c->reg.clipdistance_offset, brw_imm_d(clipdist0_offset - 6*sizeof(float)));   brw_DO(p, BRW_EXECUTE_1);   {      /* if (planemask & 1)       */      brw_AND(p, v1_null_ud, c->reg.planemask, brw_imm_ud(1));      brw_inst_set_cond_modifier(p->devinfo, brw_last_inst, BRW_CONDITIONAL_NZ);      brw_IF(p, BRW_EXECUTE_1);      {         brw_AND(p, v1_null_ud, c->reg.vertex_src_mask, brw_imm_ud(1));         brw_inst_set_cond_modifier(p->devinfo, brw_last_inst, BRW_CONDITIONAL_NZ);         brw_IF(p, BRW_EXECUTE_1);         {            /* user clip distance: just fetch the correct float from each vertex */            struct brw_indirect temp_ptr = brw_indirect(7, 0);            brw_ADD(p, get_addr_reg(temp_ptr), get_addr_reg(vtx0), c->reg.clipdistance_offset);            brw_MOV(p, c->reg.dp0, deref_1f(temp_ptr, 0));            brw_ADD(p, get_addr_reg(temp_ptr), get_addr_reg(vtx1), c->reg.clipdistance_offset);            brw_MOV(p, c->reg.dp1, deref_1f(temp_ptr, 0));         }         brw_ELSE(p);         {            /* fixed plane: fetch the hpos, dp4 against the plane. */            if (c->key.nr_userclip)               brw_MOV(p, c->reg.plane_equation, deref_4f(plane_ptr, 0));            else               brw_MOV(p, c->reg.plane_equation, deref_4b(plane_ptr, 0));            brw_DP4(p, vec4(c->reg.dp0), deref_4f(vtx0, hpos_offset), c->reg.plane_equation);            brw_DP4(p, vec4(c->reg.dp1), deref_4f(vtx1, hpos_offset), c->reg.plane_equation);         }//.........这里部分代码省略.........

void nuiDrawContext::DrawGradient(const nuiGradient& rGradient, const nuiRect& rEnclosingRect, nuiSize x1, nuiSize y1, nuiSize x2, nuiSize y2){  nuiVector2 vec(x2 - x1, y2 - y1);  nuiVector2 para(-vec[1], vec[0]);  nuiVector2 vec1(vec);  nuiVector2 para1(para);  vec1.Normalize();  para1.Normalize();  // What Quadrant are we in?:  //         |  //     a   |   b  //         |  //  ----------------  //         |  //     c   |   d  //         |  float xa, xb, xc, xd;  float ya, yb, yc, yd;  float x, y;  float xp, yp;  float xx, yy;  float xxp, yyp;  xa = xc = rEnclosingRect.Left();  xb = xd = rEnclosingRect.Right();  ya = yb = rEnclosingRect.Top();  yc = yd = rEnclosingRect.Bottom();  if (x1 < x2)  {    // Go from a to d or c to b    if (y1 == y2)    {      x  = xa; y  = ya;      xp = xc; yp = yc;      xx = xd; yy = yd;      xxp= xb; yyp= yb;    }    else if (y1 < y2)    {      // a to d      IntersectLines(xa,ya, para1[0], para1[1], xb, yb, vec1[0], vec1[1], x, y);      IntersectLines(xa,ya, para1[0], para1[1], xc, yc, vec1[0], vec1[1], xp, yp);      IntersectLines(xd,yd, para1[0], para1[1], xc, yc, vec1[0], vec1[1], xx, yy);      IntersectLines(xd,yd, para1[0], para1[1], xb, yb, vec1[0], vec1[1], xxp, yyp);    }    else    {      // c to d      IntersectLines(xc,yc, para1[0], para1[1], xa, ya, vec1[0], vec1[1], x, y);      IntersectLines(xc,yc, para1[0], para1[1], xd, yd, vec1[0], vec1[1], xp, yp);      IntersectLines(xb,yb, para1[0], para1[1], xd, yd, vec1[0], vec1[1], xx, yy);      IntersectLines(xb,yb, para1[0], para1[1], xa, ya, vec1[0], vec1[1], xxp, yyp);    }  }  else  {    if (y1 == y2)    {      x  = xd; y  = yd;      xp = xb; yp = yb;      xx = xa; yy = ya;      xxp= xc; yyp= yc;    }    else if (y1 < y2)    {      // b to c      IntersectLines(xb,yb, para1[0], para1[1], xd, yd, vec1[0], vec1[1], x, y);      IntersectLines(xb,yb, para1[0], para1[1], xa, ya, vec1[0], vec1[1], xp, yp);      IntersectLines(xc,yc, para1[0], para1[1], xa, ya, vec1[0], vec1[1], xx, yy);      IntersectLines(xc,yc, para1[0], para1[1], xd, yd, vec1[0], vec1[1], xxp, yyp);    }    else    {      // d to a      IntersectLines(xd,yd, para1[0], para1[1], xc, yc, vec1[0], vec1[1], x, y);      IntersectLines(xd,yd, para1[0], para1[1], xb, yb, vec1[0], vec1[1], xp, yp);      IntersectLines(xa,ya, para1[0], para1[1], xb, yb, vec1[0], vec1[1], xx, yy);      IntersectLines(xa,ya, para1[0], para1[1], xc, yc, vec1[0], vec1[1], xxp, yyp);    }  }  float startx,starty;  float startxp,startyp;  float stopx,stopy;  float stopxp,stopyp;  if (y1 != y2)  {    IntersectLines(x1, y1, para1[0], para1[1], x,  y,  vec1[0], vec1[1], startx,  starty);    IntersectLines(x1, y1, para1[0], para1[1], xp, yp, vec1[0], vec1[1], startxp, startyp);    IntersectLines(x2, y2, para1[0], para1[1], x,  y,  vec1[0], vec1[1], stopx,   stopy);    IntersectLines(x2, y2, para1[0], para1[1], xp, yp, vec1[0], vec1[1], stopxp,  stopyp);  }  else  {    startx  = x1; starty  = y;    startxp = x1; startyp = yp;    stopx   = x2; stopy   = y;//.........这里部分代码省略.........

boolMAST::HeatConductionElementBase::internal_residual (bool request_jacobian,                                                    RealVectorX& f,                                                    RealMatrixX& jac) {        const std::vector<Real>& JxW           = _fe->get_JxW();    const std::vector<libMesh::Point>& xyz = _fe->get_xyz();    const unsigned int    n_phi  = _fe->n_shape_functions(),    dim    = _elem.dim();        RealMatrixX    material_mat   = RealMatrixX::Zero(dim, dim),    dmaterial_mat  = RealMatrixX::Zero(dim, dim), // for calculation of Jac when k is temp. dep.    mat_n2n2       = RealMatrixX::Zero(n_phi, n_phi);    RealVectorX    vec1     = RealVectorX::Zero(1),    vec2_n2  = RealVectorX::Zero(n_phi),    flux     = RealVectorX::Zero(dim);        std::auto_ptr<MAST::FieldFunction<RealMatrixX> > conductance =    _property.thermal_conductance_matrix(*this);        libMesh::Point p;    std::vector<MAST::FEMOperatorMatrix> dBmat(dim);    MAST::FEMOperatorMatrix Bmat; // for calculation of Jac when k is temp. dep.        for (unsigned int qp=0; qp<JxW.size(); qp++) {        // initialize the Bmat operator for this term        _initialize_mass_fem_operator(qp, Bmat);        Bmat.right_multiply(vec1, _sol);        if (_active_sol_function)            dynamic_cast<MAST::MeshFieldFunction<RealVectorX>*>            (_active_sol_function)->set_element_quadrature_point_solution(vec1);        _local_elem->global_coordinates_location(xyz[qp], p);                (*conductance)(p, _time, material_mat);        _initialize_flux_fem_operator(qp, dBmat);                // calculate the flux for each dimension and add its weighted        // component to the residual        flux.setZero();        for (unsigned int j=0; j<dim; j++) {            dBmat[j].right_multiply(vec1, _sol);        // dT_dxj                        for (unsigned int i=0; i<dim; i++)                flux(i) += vec1(0) * material_mat(i,j); // q_i = k_ij dT_dxj        }        // now add to the residual vector        for (unsigned int i=0; i<dim; i++) {            vec1(0)  = flux(i);            dBmat[i].vector_mult_transpose(vec2_n2, vec1);            f += JxW[qp] * vec2_n2;        }                if (request_jacobian) {                        // Jacobian contribution from int_omega dB_dxi^T k_ij dB_dxj            for (unsigned int i=0; i<dim; i++)                for (unsigned int j=0; j<dim; j++) {                                        dBmat[i].right_multiply_transpose(mat_n2n2, dBmat[j]);                    jac += JxW[qp] * material_mat(i,j) * mat_n2n2;                }                        // Jacobian contribution from int_omega dB_dxi dT_dxj dk_ij/dT B            if (_active_sol_function) {                // get derivative of the conductance matrix wrt temperature                conductance->derivative(MAST::PARTIAL_DERIVATIVE,                                        *_active_sol_function,                                        p,                                        _time, dmaterial_mat);                                for (unsigned int j=0; j<dim; j++) {                    dBmat[j].right_multiply(vec1, _sol);  // dT_dxj                    for (unsigned int i=0; i<dim; i++)                        if (dmaterial_mat(i,j) != 0.) { // no need to process for zero terms                            // dB_dxi^T B                            dBmat[i].right_multiply_transpose(mat_n2n2, Bmat);                            // dB_dxi^T (dT_dxj dk_ij/dT) B                            jac += JxW[qp] * vec1(0) * dmaterial_mat(i,j) * mat_n2n2;                        }                }            }        }    }        if (_active_sol_function)        dynamic_cast<MAST::MeshFieldFunction<RealVectorX>*>        (_active_sol_function)->clear_element_quadrature_point_solution();    return request_jacobian;//.........这里部分代码省略.........

boolMAST::HeatConductionElementBase::velocity_residual (bool request_jacobian,                                                    RealVectorX& f,                                                    RealMatrixX& jac_xdot,                                                    RealMatrixX& jac) {    MAST::FEMOperatorMatrix Bmat;        const std::vector<Real>& JxW                 = _fe->get_JxW();    const std::vector<libMesh::Point>& xyz       = _fe->get_xyz();        const unsigned int    n_phi      = _fe->n_shape_functions(),    dim        = _elem.dim();        RealMatrixX    material_mat    = RealMatrixX::Zero(dim, dim),    mat_n2n2        = RealMatrixX::Zero(n_phi, n_phi);    RealVectorX    vec1    = RealVectorX::Zero(1),    vec2_n2 = RealVectorX::Zero(n_phi);        std::auto_ptr<MAST::FieldFunction<RealMatrixX> > capacitance =    _property.thermal_capacitance_matrix(*this);        libMesh::Point p;        for (unsigned int qp=0; qp<JxW.size(); qp++) {                _initialize_mass_fem_operator(qp, Bmat);        Bmat.right_multiply(vec1, _sol);               //  B * T                if (_active_sol_function)            dynamic_cast<MAST::MeshFieldFunction<RealVectorX>*>            (_active_sol_function)->set_element_quadrature_point_solution(vec1);        _local_elem->global_coordinates_location(xyz[qp], p);                (*capacitance)(p, _time, material_mat);                Bmat.right_multiply(vec1, _vel);               //  B * T_dot        Bmat.vector_mult_transpose(vec2_n2, vec1);     //  B^T * B * T_dot                f      += JxW[qp] * material_mat(0,0) * vec2_n2; // (rho*cp)*JxW B^T B T_dot                if (request_jacobian) {                        Bmat.right_multiply_transpose(mat_n2n2, Bmat);  // B^T B            jac_xdot += JxW[qp] * material_mat(0,0) * mat_n2n2;  // B^T B * JxW (rho*cp)                        // Jacobian contribution from int_omega B T d(rho*cp)/dT B            if (_active_sol_function) {                // get derivative of the conductance matrix wrt temperature                capacitance->derivative(MAST::PARTIAL_DERIVATIVE,                                        *_active_sol_function,                                        p,                                        _time, material_mat);                                if (material_mat(0,0) != 0.) { // no need to process for zero terms                                        // B^T (T d(rho cp)/dT) B                    jac += JxW[qp] * vec1(0) * material_mat(0,0) * mat_n2n2;                }            }        }    }            if (_active_sol_function)        dynamic_cast<MAST::MeshFieldFunction<RealVectorX>*>        (_active_sol_function)->clear_element_quadrature_point_solution();    return request_jacobian;}

/*********************************************************************** * Output clipped polygon as an unfilled primitive: */static void emit_lines(struct brw_clip_compile *c,		       bool do_offset){   struct brw_compile *p = &c->func;   const struct brw_context *brw = p->brw;   struct brw_indirect v0 = brw_indirect(0, 0);   struct brw_indirect v1 = brw_indirect(1, 0);   struct brw_indirect v0ptr = brw_indirect(2, 0);   struct brw_indirect v1ptr = brw_indirect(3, 0);   /* Need a seperate loop for offset:    */   if (do_offset) {      brw_MOV(p, c->reg.loopcount, c->reg.nr_verts);      brw_MOV(p, get_addr_reg(v0ptr), brw_address(c->reg.inlist));      brw_DO(p, BRW_EXECUTE_1);      {	 brw_MOV(p, get_addr_reg(v0), deref_1uw(v0ptr, 0));	 brw_ADD(p, get_addr_reg(v0ptr), get_addr_reg(v0ptr), brw_imm_uw(2));		 apply_one_offset(c, v0);		 brw_ADD(p, c->reg.loopcount, c->reg.loopcount, brw_imm_d(-1));         brw_inst_set_cond_modifier(brw, brw_last_inst, BRW_CONDITIONAL_G);      }      brw_WHILE(p);      brw_inst_set_pred_control(brw, brw_last_inst, BRW_PREDICATE_NORMAL);   }   /* v1ptr = &inlist[nr_verts]    * *v1ptr = v0    */   brw_MOV(p, c->reg.loopcount, c->reg.nr_verts);   brw_MOV(p, get_addr_reg(v0ptr), brw_address(c->reg.inlist));   brw_ADD(p, get_addr_reg(v1ptr), get_addr_reg(v0ptr), retype(c->reg.nr_verts, BRW_REGISTER_TYPE_UW));   brw_ADD(p, get_addr_reg(v1ptr), get_addr_reg(v1ptr), retype(c->reg.nr_verts, BRW_REGISTER_TYPE_UW));   brw_MOV(p, deref_1uw(v1ptr, 0), deref_1uw(v0ptr, 0));   brw_DO(p, BRW_EXECUTE_1);   {      brw_MOV(p, get_addr_reg(v0), deref_1uw(v0ptr, 0));      brw_MOV(p, get_addr_reg(v1), deref_1uw(v0ptr, 2));      brw_ADD(p, get_addr_reg(v0ptr), get_addr_reg(v0ptr), brw_imm_uw(2));      /* draw edge if edgeflag != 0 */      brw_CMP(p,	      vec1(brw_null_reg()), BRW_CONDITIONAL_NZ,	      deref_1f(v0, brw_varying_to_offset(&c->vue_map,                                                 VARYING_SLOT_EDGE)),	      brw_imm_f(0));      brw_IF(p, BRW_EXECUTE_1);      {	 brw_clip_emit_vue(c, v0, BRW_URB_WRITE_ALLOCATE_COMPLETE,                           (_3DPRIM_LINESTRIP << URB_WRITE_PRIM_TYPE_SHIFT)                           | URB_WRITE_PRIM_START);	 brw_clip_emit_vue(c, v1, BRW_URB_WRITE_ALLOCATE_COMPLETE,                           (_3DPRIM_LINESTRIP << URB_WRITE_PRIM_TYPE_SHIFT)                           | URB_WRITE_PRIM_END);      }      brw_ENDIF(p);      brw_ADD(p, c->reg.loopcount, c->reg.loopcount, brw_imm_d(-1));      brw_inst_set_cond_modifier(brw, brw_last_inst, BRW_CONDITIONAL_NZ);   }   brw_WHILE(p);   brw_inst_set_pred_control(brw, brw_last_inst, BRW_PREDICATE_NORMAL);}

voidvec4_generator::generate_gs_set_channel_masks(struct brw_reg dst,        struct brw_reg src){    /* From p21 of volume 4 part 2 of the Ivy Bridge PRM (2.4.3.1 Message     * Header: M0.5):     *     *     15 Vertex 1 DATA [3] / Vertex 0 DATA[7] Channel Mask     *     *        When Swizzle Control = URB_INTERLEAVED this bit controls Vertex 1     *        DATA[3], when Swizzle Control = URB_NOSWIZZLE this bit controls     *        Vertex 0 DATA[7].  This bit is ANDed with the corresponding     *        channel enable to determine the final channel enable.  For the     *        URB_READ_OWORD & URB_READ_HWORD messages, when final channel     *        enable is 1 it indicates that Vertex 1 DATA [3] will be included     *        in the writeback message.  For the URB_WRITE_OWORD &     *        URB_WRITE_HWORD messages, when final channel enable is 1 it     *        indicates that Vertex 1 DATA [3] will be written to the surface.     *     *        0: Vertex 1 DATA [3] / Vertex 0 DATA[7] channel not included     *        1: Vertex DATA [3] / Vertex 0 DATA[7] channel included     *     *     14 Vertex 1 DATA [2] Channel Mask     *     13 Vertex 1 DATA [1] Channel Mask     *     12 Vertex 1 DATA [0] Channel Mask     *     11 Vertex 0 DATA [3] Channel Mask     *     10 Vertex 0 DATA [2] Channel Mask     *      9 Vertex 0 DATA [1] Channel Mask     *      8 Vertex 0 DATA [0] Channel Mask     *     * (This is from a section of the PRM that is agnostic to the particular     * type of shader being executed, so "Vertex 0" and "Vertex 1" refer to     * geometry shader invocations 0 and 1, respectively).  Since we have the     * enable flags for geometry shader invocation 0 in bits 3:0 of DWORD 0,     * and the enable flags for geometry shader invocation 1 in bits 7:0 of     * DWORD 4, we just need to OR them together and store the result in bits     * 15:8 of DWORD 5.     *     * It's easier to get the EU to do this if we think of the src and dst     * registers as composed of 32 bytes each; then, we want to pick up the     * contents of bytes 0 and 16 from src, OR them together, and store them in     * byte 21.     *     * We can do that by the following EU instruction:     *     *     or(1) dst.21<1>UB src<0,1,0>UB src.16<0,1,0>UB { align1 WE_all }     *     * Note: this relies on the source register having zeros in (a) bits 7:4 of     * DWORD 0 and (b) bits 3:0 of DWORD 4.  We can rely on (b) because the     * source register was prepared by GS_OPCODE_PREPARE_CHANNEL_MASKS (which     * shifts DWORD 4 left by 4 bits), and we can rely on (a) because prior to     * the execution of GS_OPCODE_PREPARE_CHANNEL_MASKS, DWORDs 0 and 4 need to     * contain valid channel mask values (which are in the range 0x0-0xf).     */    dst = retype(dst, BRW_REGISTER_TYPE_UB);    src = retype(src, BRW_REGISTER_TYPE_UB);    brw_push_insn_state(p);    brw_set_access_mode(p, BRW_ALIGN_1);    brw_set_mask_control(p, BRW_MASK_DISABLE);    brw_OR(p, suboffset(vec1(dst), 21), vec1(src), suboffset(vec1(src), 16));    brw_pop_insn_state(p);}

void BaseVH::filterContoursByEdgeAngle( std::vector< std::vector< cv::Point2f > > &contours , double angleThreshold ){  int numContours = contours.size();  double thresholdVal = - std::cos( ( angleThreshold / 180 ) * CV_PI );  std::vector< std::vector< cv::Point2f > > filteredContours;  double averageContourLength = 0;  int totalNumEdges = 0;  for( int cc = 0; cc < numContours; cc++ )  {	  int numEdges = contours[ cc ].size();	  std::vector< cv::Point2f > filteredContour;	  int id1 = -1 , id2 = -1 , id3 = -1;	  for( int ee = 0; ee < numEdges; ee++ )	  {		  cv::Point2f &pt1 = contours[ cc ][ ( ee - 1 + numEdges ) % numEdges ];		  cv::Point2f &pt2 = contours[ cc ][ ee ];		  tr::Vector2f vec1( pt2.x - pt1.x , pt2.y - pt1.y );		  averageContourLength += vec1.norm();	  }	  totalNumEdges += numEdges;  }  if( totalNumEdges < 4 )	  return;  averageContourLength /= totalNumEdges;  double lengthThreshold = 1e-4 * averageContourLength;  for( int cc = 0; cc < numContours; cc++ )  {	  int numEdges = contours[ cc ].size();	  std::vector< cv::Point2f > filteredContour;	  int id1 = -1 , id2 = -1 , id3 = -1;	  cv::Point2f prevPoint( 0 , 0 );	  for( int ee = 0; ee < numEdges; ee++ )	  {		  cv::Point2f &pt1 = contours[ cc ][ ( ee - 1 + numEdges ) % numEdges ];		  cv::Point2f &pt2 = contours[ cc ][ ee ];		  cv::Point2f &pt3 = contours[ cc ][ ( ee + 1 ) % numEdges ];		  tr::Vector2f vec1( pt2.x - pt1.x , pt2.y - pt1.y ) , vec2( pt3.x - pt2.x , pt3.y - pt2.y );		  vec1.normalize();		  vec2.normalize();  		  double val = vec1.dot( vec2 );		  if( val < thresholdVal )		  {		    continue;		  }		  if( filteredContour.size() > 0 )		  {			  tr::Vector2f vec( prevPoint.x - pt2.x , prevPoint.y - pt2.y );			  if( vec.norm() < lengthThreshold )			  {			    continue;			  }		  }		  prevPoint = pt2;		  filteredContour.push_back( pt2 );	  }	  filteredContours.push_back( filteredContour );  }  contours = filteredContours;}

/***   Since no auxiliary routines are needed, we can directly start with main().**/int main(){  typedef float       ScalarType;  /**  * Initialize OpenCL vectors:  **/  unsigned int vector_size = 10;  viennacl::vector<ScalarType>  vec1(vector_size);  viennacl::vector<ScalarType>  vec2(vector_size);  viennacl::vector<ScalarType>  result_mul(vector_size);  viennacl::vector<ScalarType>  result_div(vector_size);  /**  * Fill the operands vec1 and vec2 with some numbers.  **/  for (unsigned int i=0; i<vector_size; ++i)  {    vec1[i] = static_cast<ScalarType>(i);    vec2[i] = static_cast<ScalarType>(vector_size-i);  }  /**  * Set up the OpenCL program given in my_compute_kernel:  * A program is one compilation unit and can hold many different compute kernels.  **/  viennacl::ocl::program & my_prog = viennacl::ocl::current_context().add_program(my_compute_program, "my_compute_program");  // Note: Releases older than ViennaCL 1.5.0 required calls to add_kernel(). This is no longer needed, the respective interface has been removed.  /**  * Now we can get the kernels from the program 'my_program'.  * (Note that first all kernels need to be registered via add_kernel() before get_kernel() can be called,  * otherwise existing references might be invalidated)  **/  viennacl::ocl::kernel & my_kernel_mul = my_prog.get_kernel("elementwise_prod");  viennacl::ocl::kernel & my_kernel_div = my_prog.get_kernel("elementwise_div");  /**  * Launch the kernel with 'vector_size' threads in one work group  * Note that std::size_t might differ between host and device. Thus, a cast to cl_uint is necessary for the forth argument.  **/  viennacl::ocl::enqueue(my_kernel_mul(vec1, vec2, result_mul, static_cast<cl_uint>(vec1.size())));  viennacl::ocl::enqueue(my_kernel_div(vec1, vec2, result_div, static_cast<cl_uint>(vec1.size())));  /**  * Print the result:  **/  std::cout << "        vec1: " << vec1 << std::endl;  std::cout << "        vec2: " << vec2 << std::endl;  std::cout << "vec1 .* vec2: " << result_mul << std::endl;  std::cout << "vec1 /* vec2: " << result_div << std::endl;  std::cout << "norm_2(vec1 .* vec2): " << viennacl::linalg::norm_2(result_mul) << std::endl;  std::cout << "norm_2(vec1 /* vec2): " << viennacl::linalg::norm_2(result_div) << std::endl;  /**  *  We are already done. We only needed a few lines of code by letting ViennaCL deal with the details :-)  **/  std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;  return EXIT_SUCCESS;}

/* Use mesa's clipping algorithms, translated to GEN4 assembly. */void brw_clip_tri( struct brw_clip_compile *c ){   struct brw_compile *p = &c->func;   struct brw_indirect vtx = brw_indirect(0, 0);   struct brw_indirect vtxPrev = brw_indirect(1, 0);   struct brw_indirect vtxOut = brw_indirect(2, 0);   struct brw_indirect plane_ptr = brw_indirect(3, 0);   struct brw_indirect inlist_ptr = brw_indirect(4, 0);   struct brw_indirect outlist_ptr = brw_indirect(5, 0);   struct brw_indirect freelist_ptr = brw_indirect(6, 0);   struct brw_instruction *plane_loop;   struct brw_instruction *plane_active;   struct brw_instruction *vertex_loop;   struct brw_instruction *next_test;   struct brw_instruction *prev_test;      brw_MOV(p, get_addr_reg(vtxPrev),     brw_address(c->reg.vertex[2]) );   brw_MOV(p, get_addr_reg(plane_ptr),   brw_clip_plane0_address(c));   brw_MOV(p, get_addr_reg(inlist_ptr),  brw_address(c->reg.inlist));   brw_MOV(p, get_addr_reg(outlist_ptr), brw_address(c->reg.outlist));   brw_MOV(p, get_addr_reg(freelist_ptr), brw_address(c->reg.vertex[3]) );   plane_loop = brw_DO(p, BRW_EXECUTE_1);   {      /* if (planemask & 1)       */      brw_set_conditionalmod(p, BRW_CONDITIONAL_NZ);      brw_AND(p, vec1(brw_null_reg()), c->reg.planemask, brw_imm_ud(1));            plane_active = brw_IF(p, BRW_EXECUTE_1);      {	 /* vtxOut = freelist_ptr++ 	  */	 brw_MOV(p, get_addr_reg(vtxOut),       get_addr_reg(freelist_ptr) );	 brw_ADD(p, get_addr_reg(freelist_ptr), get_addr_reg(freelist_ptr), brw_imm_uw(c->nr_regs * REG_SIZE));	 if (c->key.nr_userclip)	    brw_MOV(p, c->reg.plane_equation, deref_4f(plane_ptr, 0));	 else	    brw_MOV(p, c->reg.plane_equation, deref_4b(plane_ptr, 0));	    	 brw_MOV(p, c->reg.loopcount, c->reg.nr_verts);	 brw_MOV(p, c->reg.nr_verts, brw_imm_ud(0));	 vertex_loop = brw_DO(p, BRW_EXECUTE_1);	 {	    /* vtx = *input_ptr;	     */	    brw_MOV(p, get_addr_reg(vtx), deref_1uw(inlist_ptr, 0));	    /* IS_NEGATIVE(prev) */	    brw_set_conditionalmod(p, BRW_CONDITIONAL_L);	    brw_DP4(p, vec4(c->reg.dpPrev), deref_4f(vtxPrev, c->offset_hpos), c->reg.plane_equation);	    prev_test = brw_IF(p, BRW_EXECUTE_1);	    {	       /* IS_POSITIVE(next)		*/	       brw_set_conditionalmod(p, BRW_CONDITIONAL_GE);	       brw_DP4(p, vec4(c->reg.dp), deref_4f(vtx, c->offset_hpos), c->reg.plane_equation);	       next_test = brw_IF(p, BRW_EXECUTE_1);	       {		  /* Coming back in.		   */		  brw_ADD(p, c->reg.t, c->reg.dpPrev, negate(c->reg.dp));		  brw_math_invert(p, c->reg.t, c->reg.t);		  brw_MUL(p, c->reg.t, c->reg.t, c->reg.dpPrev);		  /* If (vtxOut == 0) vtxOut = vtxPrev		   */		  brw_CMP(p, vec1(brw_null_reg()), BRW_CONDITIONAL_EQ, get_addr_reg(vtxOut), brw_imm_uw(0) );		  brw_MOV(p, get_addr_reg(vtxOut), get_addr_reg(vtxPrev) );		  brw_set_predicate_control(p, BRW_PREDICATE_NONE);		  brw_clip_interp_vertex(c, vtxOut, vtxPrev, vtx, c->reg.t, GL_FALSE);		  /* *outlist_ptr++ = vtxOut;		   * nr_verts++; 		   * vtxOut = 0;		   */		  brw_MOV(p, deref_1uw(outlist_ptr, 0), get_addr_reg(vtxOut));		  brw_ADD(p, get_addr_reg(outlist_ptr), get_addr_reg(outlist_ptr), brw_imm_uw(sizeof(short)));		  brw_ADD(p, c->reg.nr_verts, c->reg.nr_verts, brw_imm_ud(1));		  brw_MOV(p, get_addr_reg(vtxOut), brw_imm_uw(0) );	       }	       brw_ENDIF(p, next_test);	       	    }	    prev_test = brw_ELSE(p, prev_test);	    {	       /* *outlist_ptr++ = vtxPrev;		* nr_verts++;		*/	       brw_MOV(p, deref_1uw(outlist_ptr, 0), get_addr_reg(vtxPrev));	       brw_ADD(p, get_addr_reg(outlist_ptr), get_addr_reg(outlist_ptr), brw_imm_uw(sizeof(short)));	       brw_ADD(p, c->reg.nr_verts, c->reg.nr_verts, brw_imm_ud(1));//.........这里部分代码省略.........

Matrix<double> BspCurvBasisFuncSet::CreateMatrixIntegral(int lev, double x1, double x2) const{		KnotSet kset = KnotSet(*kts,ord,num).CreateKnotSetDeriv(lev);		Matrix<double> mat(kset.GetNum()-(ord-lev),kset.GetNum()-(ord-lev));		BspCurvBasisFuncSet basis(kset.GetKnots(),ord-lev,kset.GetNum());	// find the first basis function for which the last distinct	// knot is greater than x1	int ind1=-1;	do {		ind1++;	} while ((*(basis.b))[ind1].GetKnots()[ord-lev] <= x1);	// find the last basis function for which the first distinct	// knot is less than x2	int ind2=-1;	do {		ind2++;	} while (ind2 < kset.GetNum()-ord+lev && (*(basis.b))[ind2].GetKnots()[0] < x2);	Matrix<double> mat1(kset.GetNum()-ord+lev,kset.GetNum()-ord+lev,0.0);	for (int i=ind1; i<=ind2-1; i++)		for (int j=ind1; j<=i; j++) {								// create the two std::sets representing the two knot std::sets			Vector<double> temp1((*(basis.b))[i].GetKnots());			Vector<double> temp2((*(basis.b))[j].GetKnots());			std::set<double> s1(temp1.begin(),temp1.end());			std::set<double> s2(temp2.begin(),temp2.end());			if (*(--s2.end()) > *(s1.begin())) {				// form the intersection				std::set<double> s3;				std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(),std::inserter(s3,s3.begin()));							// if there is an intersection				if (s3.size() > 1) {					Vector<double> v(s3.size());					std::set<double>::iterator s = s3.begin();					// copy the elements into a vector					for (unsigned int k=0; k<s3.size(); k++) v[k] = *s++;									// create the compbezcurvs					Vector<BezCurv<double> > vec1(s3.size()-1);					Vector<BezCurv<double> > vec2(s3.size()-1);					BspCurv<double> b1((*(basis.b))[i].GetBspCurv()), b2((*(basis.b))[j].GetBspCurv());									// find the segments of intersection					for (unsigned int k=0; k<s3.size()-1; k++) {						int segb1 = b1.GetKnotSet().Find_segment(v[k]);						int segb2 = b2.GetKnotSet().Find_segment(v[k]);												vec1[k] = b1.GetSegment(segb1);						vec2[k] = b2.GetSegment(segb2);					}									CompBezCurv<double> cb1(vec1,s3.size()-1,v);					CompBezCurv<double> cb2(vec2,s3.size()-1,v);					CompBezCurv<double> prod = cb1.Product(cb2);									mat1[i][j] = prod.ConvertBspCurv().Integrate(x1,x2);				}			}		}		for (int i=ind1; i<=ind2-2; i++)		for (int j=i+1; j<=ind2-1; j++) mat1[i][j] = mat1[j][i];		return mat1;}

void Ellipsoid::subDraw(){	computeAssistVar();	getPrimitiveSetList().clear();	osg::ref_ptr<osg::Vec3Array> vertexArr = new osg::Vec3Array;	osg::ref_ptr<osg::Vec3Array> normalArr = new osg::Vec3Array;	setVertexArray(vertexArr);	setNormalArray(normalArr, osg::Array::BIND_PER_VERTEX);	osg::ref_ptr<osg::Vec4Array> colArr = new osg::Vec4Array();	colArr->push_back(m_color);	setColorArray(colArr, osg::Array::BIND_OVERALL);	bool isFull = osg::equivalent(m_angle, 2 * M_PI, GetEpsilon());	if (isFull)	{		m_angle = 2 * M_PI;	}	osg::Vec3 bottomNormal = -m_aLen;	bottomNormal.normalize();	osg::Quat localToWold;	localToWold.makeRotate(osg::Z_AXIS, -bottomNormal);	osg::Vec3 xVec = localToWold * osg::X_AXIS;	osg::Vec3 yVec = xVec ^ bottomNormal;	int hCount = m_bDivision;	double hIncAng = 2 * M_PI / hCount;	osg::Quat hQuat(hIncAng, bottomNormal);	int vCount = (int)ceil(m_angle / (2 * M_PI / m_aDivision));	if (vCount & 1) // 如果是奇数，则变成偶数		++vCount;	double vIncAng = m_angle / vCount;	double currAngle = m_angle / 2.0;	double b = m_bRadius;	double a = m_aLen.length();	osg::Vec3 vec1(b * sin(currAngle), 0, a * cos(currAngle));	vec1 = localToWold * vec1;	osg::Vec3 normal1(sin(currAngle) / b, 0, cos(currAngle) / a);	normal1 = localToWold * normal1;	normal1.normalize();	currAngle -= vIncAng;	osg::Vec3 vec2(b * sin(currAngle), 0, a * cos(currAngle));	vec2 = localToWold * vec2;	osg::Vec3 normal2(sin(currAngle) / b, 0, cos(currAngle) / a);	normal2 = localToWold * normal2;	normal2.normalize();	const GLint first = vertexArr->size();	for (int i = 0; i < vCount / 2; ++i)	{		osg::Vec3 hVec1 = vec1;		osg::Vec3 hVec2 = vec2;		osg::Vec3 hNormal1 = normal1;		osg::Vec3 hNormal2 = normal2;		const size_t hFirst = vertexArr->size();		for (int j = 0; j < hCount; ++j)		{			vertexArr->push_back(m_center + hVec1);			vertexArr->push_back(m_center + hVec2);			normalArr->push_back(hNormal1);			normalArr->push_back(hNormal2);			hVec1 = hQuat * hVec1;			hVec2 = hQuat * hVec2;			hNormal1 = hQuat * hNormal1;			hNormal2 = hQuat * hNormal2;		}		vertexArr->push_back((*vertexArr)[hFirst]);		vertexArr->push_back((*vertexArr)[hFirst + 1]);		normalArr->push_back((*normalArr)[hFirst]);		normalArr->push_back((*normalArr)[hFirst + 1]);		vec1 = vec2;		currAngle -= vIncAng;		vec2.set(b * sin(currAngle), 0, a * cos(currAngle));		vec2 = localToWold * vec2;		normal1 = normal2;		normal2.set(sin(currAngle) / b, 0, cos(currAngle) / a);		normal2 = localToWold * normal2;		normal2.normalize();	}	addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::QUAD_STRIP, first, vertexArr->size() - first));	if (!isFull && m_bottomVis)	{		const GLint first = vertexArr->size();		currAngle = m_angle / 2.0;		osg::Vec3 vec1(b * sin(currAngle), 0, a * cos(currAngle));		vec1 = localToWold * vec1;		osg::Vec3 vec2(0, 0, a * cos(currAngle));		vec2 = localToWold * vec2;		vertexArr->push_back(m_center + vec2);		normalArr->push_back(bottomNormal);		for (int i = 0; i < hCount; ++i)//.........这里部分代码省略.........